Guillaume Charpiat Pulsar
email: guillaume.charpiat@sophia.inria.fr

Learning Shape Metrics based on Deformations and Transport

Shape evolutions, as well as shape matchings or image segmentation with shape prior, involve the preliminary choice of a suitable metric in the space of shapes. Instead of choosing a particular one, we propose a framework to learn shape metrics from a set of examples of shapes, designed to be able to handle sparse sets of highly varying shapes, since typical shape datasets, like human silhouettes, are intrinsically high-dimensional and non-dense. We formulate the task of finding the optimal metrics on an empirical manifold of shapes as a classical minimization problem ensuring smoothness, and compute its global optimum fast.

First, we design a criterion to compute point-to-point matching between shapes which deals with topological changes. Then, given a training set of shapes, we use these matchings to transport deformations observed on any shape to any other one. Finally, we estimate the metric in the tangent space of any shape, based on transported deformations, weighted by their reliability. Experiments on difficult sets are shown, and applications are proposed.

Introduction

The notion of shape is important in many fields of computer vision, from tracking to scene understanding. As for usual object features, it can be used as a prior, as in image segmentation, or as a source of information, as in gesture classification. When image classification or segmentation tasks require high discriminative power or precision, the shape of objects naturally appear relevant to our human minds. However, shape is a complex notion which cannot be dealt with directly like a simple parameter in R n . Modeling shape manually is tedious, and one arising question is the one of learning shapes automatically.

Various statistical models of shapes exist in the literature. Most of them consist in estimating a mean pattern and characteristic deformations of a given set of shapes, under the assumption that the shape variability in the training set is not too high, in order to be able to consider sensible deformations from one shape to another one, and to compute sensible statistics on these deformations (usually with principal component analysis). Deformations can be given by the user, fully or with a few landmarks [START_REF] Taylor | Active shape models -'smart snakes[END_REF], or can be computed via a search for best diffeomorphisms [START_REF] Dryden | Statistical Shape Analysis[END_REF][START_REF] Vaillant | Statistics on diffeomorphisms via tangent space representation[END_REF] or matchings [START_REF] Gdalyahu | Flexible syntactic matching of curves and its application to automatic hierarchical classification of silhouettes[END_REF][START_REF] Ferrari | Accurate object detection with deformable shape models learnt from images[END_REF] for hand-designed metrics or criteria. They can also be computed by algorithms suited to particular shape representations [START_REF] Leventon | Statistical Shape Influence in Geodesic Active Contours[END_REF][START_REF] Cremers | Efficient kernel density estimation of shape and intensity priors for level set segmentation[END_REF] or just be distance gradients [START_REF] Charpiat | Approximations of shape metrics and application to shape warping and empirical shape statistics[END_REF]. In all cases, the statistics computed on deformations can be turned into a new metric in the tangent space of the mean pattern, acting as a deformation prior on one particular shape.

When the shape variability is too high, this paradigm fails because the automatic computation of deformations between very different shapes is not reliable, and because linear approximations of the space of shapes are not meaningful anymore. Distance-based algorithms, such as kernel methods, were proposed [START_REF] Rathi | Comparative analysis of kernel methods for statistical shape learning[END_REF][START_REF] Etyngier | Shape priors using manifold learning techniques[END_REF] to handle high variability, but at the price of considering only distances between shapes, instead of deformations, thus losing the crucial information they carry. These methods consider training sets as graphs, whose nodes are shapes and whose edges are distances (for a particular metric chosen). They assume the neighborhood of any shape to be representative of the intrinsic dimensionality of the space of shapes, and require consequently relatively high sampling densities, which are not affordable in the case of datasets with high intrinsic dimension, like human silhouettes with at least 30 degrees of freedom. Some other interesting methods are based on shape patches or parts [START_REF] Kumar | Obj cut[END_REF][START_REF] Borenstein | Class-specific, top-down segmentation[END_REF], but they sacrifice the notion of continuous, global shape. Another approach consists in searching the training set for the closest shapes to the one of interest, based on shape distances [START_REF] Rathi | A generic framework for tracking using particle filter with dynamic shape prior[END_REF] or, for videos, on time order [START_REF] Soatto | Deformotion -deforming motion, shape average and the joint registration and segmentation of images[END_REF], and then in applying classical approaches to this neighborhood only. Again, the neighborhood representativeness issue arises. Moreover, the local metrics thus computed are not guaranteed to be globally coherent as a function of the shape of interest.

This paper aims at extending the approaches based on continuous deformations, to the case of high shape variability, where the notion of mean shape is not relevant anymore, or where deformations cannot be estimated throughout the whole training set even if they still can be computed between close enough samples. In our approach, we compute point-to-point matchings between close shapes, and use them as a way to propagate information within the training set. Thus, we can transport a deformation of a shape to any other shape, with an associated reliability weight. Metrics, i.e. inner products on deformation spaces, are then estimated on the tangent space of any shape, while taking into account deformations observed at other locations, thus decreasing dramatically the sample density required and ensuring the global coherence of the manifold.

The paper is organized as follows: First, we propose a shape matching algorithm which handles topological changes. Second, we use pairwise matchings to define transport and to build a manifold-like structure. Third, we transport deformations and propose a method to estimate metrics. We show results on video datasets and, finally, in a theoretical section, we study different tracks for metric estimation, and prove that the method presented in the previous part computes the optimal metrics for a natural criterion.

Shape matching

In order to compare quantities defined on different shapes, like deformations, we need a way to transport them from shape to shape, and to do this we need point-to-point correspondences between close shapes. Since in typical shape datasets, like walking human silhouettes, topological changes are very frequent (fig. 1, left), we need a matching algorithm able to consider pairs of shapes with different topologies. The only one we found in the literature relies on successive bipartite graph matchings and spline estimations [START_REF] Belongie | Shape matching and object recognition using shape contexts[END_REF]. It imposes similar sampling rates on both shapes, which limits the precision and the regularity of the matching (fig. 1, right). Better regularity and precision (for a given time cost) will be ensured here by oversampling the target.

In the case of contours in images, a shape is a union of 1D curves. If shapes have only one connected component, i.e. if they are topologically equivalent to a circle or a segment, then dynamic time warping helps find quickly the global optima of simple matching energies [START_REF] Geiger | Dynamic programming for detecting, tracking, and matching deformable contours[END_REF][START_REF] Trouvé | Diffeomorphic matching problems in one dimension: Designing and minimizing matching functionals[END_REF][START_REF] Schmidt | Fast matching of planar shapes in sub-cubic runtime[END_REF]. In recent works, graph-cuts or minimum cycle search in graphs have also been used [START_REF] Schoenemann | Matching non-rigidly deformable shapes across images: A globally optimal solution[END_REF]. We will adapt here dynamic time warping to minimize an energy which favors smooth matchings and deals with several connected components, topological changes, and optionally with vanishing parts. Such an algorithm is needed in the sequel, but our work is not specific to the particular matching algorithm presented here.

Matching criterion

Let us start with the case of simple, closed curves A and B, seen as functions from the circle S 1 to R 2 , parameterized by their arc length s. We search for the best matching from A to B, i.e. for the best function m from S 1 to S 1 so that A ≃ B • m. Let us note f = B • m -A so that f (s) stands for the vector --------→ A (s) B (m(s)) linking a point of A to its correspondent on B. The deformation f should be as small and as smooth as possible (see figure 1 for explanation), (Middle) Matching energy explanation. We search for small linking vectors f (s), and for small difference between the two green vectors, or between the red ones, to ensure spatial coherence. These two differences have same norm: ∂f /∂s . (Right) Imposing similar sampling rates limits precision and regularity. so we minimize the following criterion over deformations f that can be written as f = B • m -A for some m:

f 2 H 1 α := S1 f (s) 2 + α ∂f ∂s 2 ds.
When shapes have thin parts, more information is required in order to distinguish nearby parallel sides. We include in the criterion the relative angle between outgoing normals at corresponding points ∠(n A(s) , n B•m(s) ) :

E match (m) = B • m -A 2 H 1 α + γ ∠(n A , n B•m ) 2 L 2
Note that we require the deformation f to be smooth, and not the parameterization correspondence m itself. Now if B as several connected components, say B = ∪ i B i , then we can still search for an optimal matching m between A and B, with m : S 1 → i S 1 possibly pointing to any of the parameterization supports of connected components B i . The matching cannot be guaranteed to be one-to-one anymore, since some parts of B may have no antecedent through m. In order to reciprocally account for points on A that cannot be matched to B (disappearing parts), we may allow a match to nothing, i.e. m(s) may have the value ∅.

Optimization with dynamic time warping

In practice, shapes are unions of polygons, and the energy can be discretized accordingly. In the case where A has only one connected component, E match can be minimized efficiently. For each vertex A(s) of A, we define the set N (s) = {∅} ∪ {s ′ s.t. B(s ′ ) -A(s) d max } of its possible matches, i.e. the set of possible values of m(s), as points of B relatively close and ∅. The consideration of a maximum distance d max speeds up the process significantly. Choosing any initial point A(s 0 = 1), and considering the set of ordered vertices along A, the problem reduces to finding an optimal function from {1, . . . , #A} to N (1) × N (2) • • • × N (#A), which can be seen as a search for an optimal path in a graph, with costs on graph edges (s, m(s)), (s + 1, m(s + 1)) derived from the energy E match . A constant high cost is assigned to edges involving two ∅, and an even higher one if involving only one ∅. This problem is solved by dynamic time warping.

If A has several components, each of them is treated independently. This process is not symmetric in the sense that the matching obtained from A to B can differ from the one from B to A. The quality and the precision of the results increase when the discretization of the target B is finer than the template A, so we oversample targets. This is corroborated by a convergence study when discretizations get finer, included in the supplementary materials.

The computational cost is low, only a fraction of a second to match shapes with hundreds of points on a standard PC. The value of E match (m) reflects the quality of the matching computed: the lower the energy is, the more similar the two shapes are, and the more reliable the matching found is. We noticed that allowing matchings to ∅ gives more accurate correspondence fields, but unluckily also less significant values E match (m) (because ∅ induces a saturation cost). Because of the importance of these reliability values in the sequel, we remove the possibility of matching to ∅.

Transport and information propagation

We now use the matching tool to define transport in training sets of shapes, and to propagate information with reliability weights. Points on shapes will now be confused with their parameterizations, so that functions can be defined on shapes rather than on the parameterizations thereof.

Local transport

Let A and B be two shapes, and m A B the matching from A to B (so that A ≃ B • m A B ). Any function h defined along B, with values in any space X , can be transported to A, or more exactly to the points of A linked with points of B, since m A B may be not one-to-one. More generally, all quantities in the sequel will be computed over matching domains. The local transport T L B A is defined by:

∀ h : B → X , T L B A (h) : A → X T L B A (h) (s) = h (m A B (s))
We can compare any two functions h A , h B defined on different shapes, by

h A -T L B A (h B ) or T L A B (h A ) -h B .
When the functions to be transported are deformations, other transports may be defined. For example, in the case of rotations, one may prefer the angle between the normal to the shape n B(s) and the vector h(s) to be kept constant during transport. However it is not obvious whether such a transport would be sensible for all deformations h and all pairs (A, B). In the sequel we keep the former transport.

Global transport

Given a training set of shapes S = (S i ), we compute all possible pairwise matchings m i j and associate to each of them the matching cost C m ij = E match (m i j ). A low cost C m ij means that the shapes S i and S j were close and that the matching m i j is reliable, whereas a high cost reveals an unsatisfying matching. For any pair (i 0 , j 0 ) we search for the best path from S i0 to S j0 in the graph whose nodes are shapes and whose edges are matching costs C m ij . We then denote by C G i0j0 the cost of this path (i 0 , i 1 , . . . , i k = j 0 ) and by T G i0 j0 the composition of local transports along it:

T G i0 j0 = T L i k-1 j0 • • • • • T L i1 i2 • T L i0 i1 .
This gives the optimal transport from S i0 to S j0 . Since E match is a quadratic energy, the optimal path will prefer series of small, reliable steps to big, uncertain jumps.

The computation of all best paths is affordable with standard shortest path algorithms. In the sequel, T i j will stand for the global transport T G i j , and w G ij = e -α T C G ij for the associated transport reliability, for a fixed positive α T . We will also denote by w L ij = e -α T C m ij the confidence in the direct matching between shapes S i and S j .

Individualized transport

One could also consider individualized transports, in the sense that the best path may depend on the transported quantity. Indeed the cost of a local transport T L i j is a sum over vertices of S i : E match (m i j ) = Si E v (m i j )(s)ds so it would make sense to consider the following individual cost, for any deformation h to be transported:

C m,ind. ij (h) = Si E v (m i j )(s) h(s) ds
so that bad matchings along a shape are not significant if they occur where there is no information to transmit.

Propagating information

This structure (T i j , w G ij ) is useful to propagate information along the training set. A local matching m i j between two close shapes can be seen as a deformation from S i to S j . Because training sets are relatively small and reliable deformations are scarce, it makes sense to complete the set of observed deformations on one shape S k by deformations observed at other locations S i for which the transport T i k is reliable. In the ideal case of a rigid object with d articulations and perfect local matchings, only one observation of each articulation moving, at any position, is sufficient to realize the full complexity of all possible movements at all positions, by transports and linear combinations.

Usual articulated human models have about d = 30 degrees of freedom, and consequently the intrinsic dimension of typical shape datasets is potentially high. This implies that there is no hope in obtaining a dense training set, even with a loose grid (say N bins for each degree of freedom, which makes N d bins), even with billions of examples. Consequently, methods involving only distances [START_REF] Rathi | Comparative analysis of kernel methods for statistical shape learning[END_REF] or nearest neighbors [START_REF] Etyngier | Shape priors using manifold learning techniques[END_REF][START_REF] Rathi | A generic framework for tracking using particle filter with dynamic shape prior[END_REF][START_REF] Soatto | Deformotion -deforming motion, shape average and the joint registration and segmentation of images[END_REF] are not likely to be successful for high d, whereas our approach based on transport of deformations reduces the number of required samples from N d to N d in the ideal case.

Suitability for video analysis

In video sequences, each frame is relatively similar to the next one, so that any temporal succession of observed shapes is a natural good path to transport through. In particular it is not possible to find a shape without close neighbors, and consequently information can be shared.

We present an example for a sequence of walking human silhouettes, from the ViHASi dataset 1 . We compute automatically pairwise matchings, and from them the best transport paths (without temporal order information). The paths obtained are strongly correlated to temporal ordering (forwards or backwards). We pick one shape, add colors randomly along it, and propagate color (as a function with values (r, g, b) in R 3 ) through transports to all training shapes, in order to visualize correspondences. The result, in figure 2 (see the supplementary materials for the whole video), is correct, validating the method. Of course a few matching mistakes are sometimes observed when auto-occultation happens, but these errors are few, and they are not propagated because these mismatchings cost more. Considering individualized transports (part 2.3) would allow the transmission of information not related to the precise location where difficulties occur. For instance a hand gesture could be transmitted to another silhouette with similar arm positions even if the legs have been crossed in the meanwhile.

Learning with shapes and transport

With such a structure, one can learn functions from a space of shapes S to any vector space X , in particular to spaces of functions defined on shapes. For example, let us consider a training set of shapes with appearance (S i , A i ) ∈ S × F(S i R n ). The appearance A i (s) could be, for any point s of any shape S i , an image patch centered on s taken from the image from which S i was segmented. The transport structure transmits examples of patches A j (m i j (s)) observed at similar locations on other shapes, and thus local statistics on patches can be performed, leading to a prior of the image given the shape. Thus, in a Bayesian framework, this structure allows the learning of segmentation/detection criteria with shape and appearance priors.

For scene understanding purposes, given videos where a few objects interact, one could also learn which parts of 1 http://dipersec.king.ac.uk/VIHASI/ objects interact (by propagating contact locations through the set of shapes), and how (by propagating gestures). The perspectives opened by this kind of framework are wide.

Learning metrics

Using previous sections, we now learn metrics, show how to use them as priors, and how to learn the whole structure.

Metric estimation with weighted

H 1 α PCA
Given a training set of shapes S = (S i ), we are now able to compute matchings m i j between close shapes, which can be seen as deformations

f i j = S j • m i j -S i .
We are also able to transport these deformations to any other shape S k :

f k i j = T i k (f i j
) with a weight combining reliability about deformation and transport computation (section 2.2):

w k i j = w L ij w G ik .
We now estimate the metric in the tangent space of the shape S k , i.e. we search for a relevant inner product in the space of deformations that can be applied to S k , based on the set of weighted deformations (f k i j , w k i j ). Principal Component Analysis (PCA) seems a reasonable way to compute statistics but we need to adapt it to probability weights and to the H 1 α product, which favors smoothness and is more coherent with the matching energy E match than the standard L 2 product. PCA is derived from an energy minimization problem: the search for the best orthonormal axes e n to project data. Here, the projection error to be minimized is:

inf en|e n ′ H 1 α =δ n=n ′ i,j w k i j f k i j - n f k i j |e n H 1 α e n 2 H 1 α
This is equivalent to the maximization problem:

sup en|e n ′ H 1 α =δ n=n ′ n i,j w k i j f k i j |e n 2 H 1 α
and to: sup

en|e n ′ H 1 α =δ n=n ′ n
e n HF H e n where F = i,j w k i j f k i j ⊗ f k i j is the weighted covariance matrix, and where H = Id -α∆ is the symmetric definite operator such that a |b H 1 α = H a |b L 2 . If we note d n = H 1/2 e n then the problem becomes: sup

dn|d n ′ L 2 =δ n=n ′ n d n H 1/2 F H 1/2 d n
so that the optimal d n are the eigenvectors of H 1/2 F H 1/2 with highest eigenvalues. As with usual PCA, we diagonalize the weighted correlation matrix M instead, given by:

M (i,j),(i ′ ,j ′ ) = w k i j f k i j w k i ′ j ′ f k i ′ j ′ H 1 α .
Let γ n be the eigenvectors of M , and λ n the eigenvalues.

One can prove that

d n = ij γ (i,j) n H 1/2 w k i j f k i j so that e n = ij γ (i,j) n w k i j f k i j
Thus, PCA of the set of weighted deformations (f k i j , w k i j ) in the tangent space of a shape S k leads to modes of deformation e n , with eigenvalues λ n , computed easily from the correlation matrix M . Based on the associated Mahalanobis distance, we set the natural inner product P between any two deformations f 1 and f 2 of S k :

f 1 | f 2 P = n 1 λ 2 n f 1 |e n H 1 α e n | f 2 H 1 α .
In practice we replace λ n by max(λ n , λ noise ) for a chosen level of noise. In case the distributions obtained along significant eigenmodes would not be Gaussian, we build histograms along eigenmodes, i.e. of f k i j | e n H 1 α . The probability distribution rebuilt from histograms, as if eigenmodes were independent, was found to be relatively close to the real one in many cases, but biases between eigenmodes may also be observed, especially for small training sets.

In figure 3 we consider a complex dancing silhouette sequence 2 with high variability and fast moves, and show first eigenmodes and histograms computed for a few shapes. Since linear combinations of first eigenmodes are most probable deformations, they can be seen as deformation priors. The priors obtained here are sensible, intuitively related to articulations and cloth moves, while the usual mean-andmodes model performs poorly and kernel methods do not lead to explicit deformation priors. The advantages over neighborhood-based methods are explained in figure 4.

An example of how to use the learned metric

The metric learned can be used as a prior on shape matching. Let A be a shape from the training set, to be matched, and B the new target. Any possible matching f = B •m-A is the sum of its projection p(f ) = n a n e n on modes estimated in the tangent space of shape A, and of the remaining part noise(f ). From the metric P estimated before, we can derive a prior in H 1 α , which associates to f its cost p(f

) 2 P + 1 λ 2 noise noise(f ) 2 H 1 α
. Histograms along first eigenmodes can also be used, their negative log-likelihood be turned into a cost, and a similar noise term be added. In both cases we want to minimize an energy of the form C((a n )) + 1

λ 2 noise noise(f ) 2 H 1 α ,
which can be expressed as: 2 From Grimage platform, https://charibdis.inrialpes.fr

inf m, - → a C( -→ a ) + 1 λ 2 noise B • m -A - n a n e n 2 H 1 α
Given any -→ a , the optimal m can be found by the method described in section 1. A classical gradient descent can then be performed on -→ a . A few steps only are needed, the process consists mainly in cutting the part of B • m -A which belongs to the span of the eigenmodes to add it to p(f ).

Learning the whole structure: second pass

The process to estimate metrics from a training set of shapes consists in three steps: computation of matchings between close shapes (section 1), transport of deformations (section 2), and turning statistics on deformations with reliability weights into metrics (section 3). If we re-run the whole process a second time, we can replace the shape matching algorithm of section 1 by the matching prior deriving from the learned metric (section 3.2). In the same spirit, in section 2, we could choose the geodesic transport associated to the learned metric. These geodesics would be easy to obtain, with a path-straightening method, since we already know point-to-point matchings as well as the metric. Thus, the choices made in these two sections can be seen as reasonable initializations aimed to be replaced with learned quantities in a second pass. Concerning the third section, one could wonder whether there could be other ways to estimate metrics, given matchings and transports. This is the subject of next part, where we will show that the method we presented already computes the optimal metrics.

Criteria on shape metrics and optima

Given a training set of shapes S = (S i ), and, for any shape S k , an empirical distribution D emp of transported deformations weighted by their reliability, one can wonder what the possible ways to estimate metrics are, whether there would be an objective criterion to assess how much a metric is suited to the set of shapes, and whether it is possible to find the optimal metrics.

Criterion in one tangent space

Let us consider first the case of the tangent space to one shape only. We are given a set of deformations f j in this tangent space T , with probability weights w j whose total sum is 1, i.e. we are given the empirical distribution:

D emp = j w j δ fj
where δ • are Dirac peaks, and we would like to find a suitable inner product P for T . To any inner product P can be associated a probability distribution over deformations f :

D P (f ) ∝ e -f 2 P Figure 3
. Dancing sequence (9s, 24Hz). Top & middle: first modes of deformation for various postures of the dancing sequence. Each mode is drawn twice, with amplitude ±λn, and associated histogram is shown. Note how the modes are sensible, related to articulations (arms, legs, elbows, ...) or dress moves. Full resolution images can be found in the supplementary materials. Bottom left: some frames of the video sequence; right: mean and first modes with the classical PCA approach on level-sets [START_REF] Leventon | Statistical Shape Influence in Geodesic Active Contours[END_REF][START_REF] Cremers | Efficient kernel density estimation of shape and intensity priors for level set segmentation[END_REF] : limbs are not correctly treated. up to a normalizing constant. Let us restrict P to be continuous with respect to an inner product P 0 proposed by default, for instance H 1 . The class of all such P is still huge, and, thanks to Riesz representation theorem, for any such P there exists a linear symmetric continuous operator A s.t.:

∀ f 1 , f 2 ∈ T, f 1 | f 2 P = A f 1 | f 2 P0 .
Since such an operator can be diagonalized, there exists an orthonormal basis (e n ) for P 0 , and real, positive coefficients (α n ) such that A = n α n e n ⊗ e n , and consequently:

∀ f 1 , f 2 ∈ T, f 1 | f 2 P = n α n f 1 |e n P0 e n | f 2 P0 which implies ∀ f ∈ T, f 2 P = n α n f | e n 2 P0
so that the associated distribution

D P (f ) := n α n π 1 2 e -αn f |en 2 P 0
is Gaussian. Reciprocally, any Gaussian distribution relates to a definite positive quadratic form, i.e. an inner product on T . Thus, a search over probability distributions derived from inner products is a search over Gaussian distributions. We would like the inner product P to be relevant to the set of deformations f j . One possible way is to pick the one whose associated distribution D P is the closest to D emp . Proposition 1. The inner product P which leads to the probability distribution D P the closest to the empirical distribution D emp = j w j δ fj for the Kullback-Leibler divergence, is the one obtained by weighted PCA on (f j , w j ).

Proof. The Kullback-Leibler divergence between any two probability distributions p 1 and p 2 is defined by:

KL(p 2 |p 1 ) = p 1 ln p 1 p 2 .
Minimizing the Kullback-Leibler divergence between p 1 and p 2 with respect to p 2 consequently leads to the minimization of E(p 2 |p 1 ) =p 1 ln p 2 . In our case this gives:

E(D P |D emp ) = - j w j ln D P (f j ) (1)
= j n w j α n f j | e n 2 P0 - 1 2 ln α n + 1 2 ln π .
If we denote by F the covariance matrix j w j f j ⊗ f j , the energy becomes, up to a constant:

n α n e n |F | e n P0 - 1 2 ln α n (2)
The minimization with respect to α n gives:

∂ αn E = e n |F | e n P0 - 1 2α n = 0. (3)
At the minimum, the derivative with respect to the unitnormed deformation e n is 0 in all directions except for pure norm variation:

∂ en E = 2 α n F e n ∝ e n (4)
which implies that e n is an eigenvector of F , say with eigenvalue λ n . Together with (3) it gives:

α n = 1 2λ 2 n
and consequently the optimal inner product that induces the closest distribution to D emp is the one related to the norm: This is, up to a constant factor, the Mahalanobis distance associated to weighted PCA on (f j , w j ), which is precisely the algorithm developed in section 3.1 with P 0 = H 1 α .

f 2 P = 1 2 n f | e n 2 P0 λ 2 n .
One might however wonder whether minimizing the Kullback-Leibler distance to a sum of Dirac peaks makes sense. Luckily, the previous proposition can be extended to the case of symmetric translation-invariant unit-mass kernels K(• -•) defined on the space T of deformations. We replace D emp by the kernel-smoothed empirical distribution

D K emp (f ) = j w j K(f j -f ).
Note : The family (f j ) is finite, so we work in a finitedimensioned subspace of the tangent space T , and K can be understood, in the simple case, as a real function multiplied by the usual Lebesgue measure d f . In the infinite dimension case, K cannot be isotropic (because it has finite mass).

Proposition 2. The inner product P which leads to the probability distribution D P the closest to the empirical distribution D K emp = j w j K(f j -•) for the Kullback-Leibler divergence, is obtained by diagonalization of the sum of the correlation matrix F and the second moment of K.

Proof. Full details in the supplementary materials. After computations, we find a similar expression to (2) except that F is replaced by

F + M K where M K = T f ⊗ f K(f ).
Note that when the kernel K gets closer to a Dirac peak, M K gets closer to 0, and we obtain proposition 1 again.

We have defined a criterion, based on the Kullback-Leibler divergence, to quantify how suitable an inner product is for a tangent space given with an empirical distribution of deformations, and we have shown how to compute the optimal one. Next sections face coherency issues.

No best smooth direction field

We would like to compute a suitable inner product P i for each tangent space T i as previously, but in a coherent way: we would like P i to be close to P k if shapes S i and S k are close. One approach would be to compute a joint PCA in the tangent spaces T i of all shapes S i simultaneously, with a regularity criterion imposing that, after transport, the eigenmodes e i n in T i should not differ too much from the ones e k n in T k (with a weight w G ik ). In the continuous case, this problem can be stated as searching for smooth vector fields e n over a manifold S. But the hairy ball theorem tells us that even in the simple case where the manifold is a sphere, there exists no non-vanishing continuous tangent vector field on the sphere. Which means that there are manifolds for which we cannot find smooth fields e n whose norm is never 0, so that global modes of deformation do not always exist.

Criterion for a smooth metric

In fact we do not need the eigenmodes e k n to be smooth with respect to the shape S k , we only need the probability distributions D P k related to them to be smooth. One way to ensure this consists in requiring the distribution D P k to be close not only to the empirical distribution D emp k in the tangent space of S k , but also to the transported empirical distributions T i k (D emp i ) from neighboring shapes S i , with a weight w G ik depending on transport reliability. Proposition 3. The metric computed in part 3.1 is the optimal metric deriving from an inner product, for the criterion:

i,k w G ik KL D P k T i k (D emp i )
where

D emp i = j w L ij δ fi j and T i k (δ f ) = δ T i k (f ) . Proof. This criterion rewrites ijk w L ij w G ik ln D P k (f k i j ) which is k KL(D P k |D T emp k )
, a sum of independent terms, where D T emp k = i,j w k i j δ f k i j is the empirical distribution of transported deformations considered in section 3.1. The optimal P k are given by proposition 1 applied independently to each tangent space T k with the distribution D T emp k . Which is precisely the content of section 3.1.

Criterion for smooth probability distributions

We could also ask for smooth distributions with an explicit regularizer term. For the sake of simplicity, let us assume that the tangent spaces T i are finite-dimensioned, so that probability distributions D Pi are just functions g i defined over T i (times the Lebesgue measure). Similarly we denote by g 0 i the density functions related to the empirical distributions D K emp i (smoothed by a kernel in order to avoid Dirac peaks). Let us consider the usual L 2 norm between these density functions over T i , and denote by T i j any choice of transport of functions defined over T i , to T j . A natural criterion to minimize would be:

E ′ (g) = i g i -g 0 i 2 L 2 (Ti) + ij w ij T i j (g i )-g j 2 L 2 (Tj )
so that the desired distributions g i are close to the empirical ones g 0 i observed in the same tangent space T i , but also so that they do not vary much when transported to close shapes S j . At the minimum of E ′ , we have: ∀i, ∂ gi E ′ = 0 = g ig 0 i + j w ij (T i j (g i )g j )T * i j + w ji (g i -T j i (g j ))

where T * i j is the adjoint of the linear application T i j . This linear system in g can be rewritten as Ag = g 0 , where A is a matrix of linear operators:

A ii = 1 + j w ij T * i j T i j + w ji A ij = -w ij T * i jw ji T j i for i = j

In fact A = Id + ε∆ where ∆ is the usual graph Laplacian, but with transports since one cannot compare directly quantities defined on different tangent spaces, and ε is related to the norm of w. Thus, A is symmetric positive definite and g = A -1 g 0 = (Id +ε∆) -1 g 0 ≃ (Id -ε∆)g 0 ≃ N ε * g 0 .

Thus, the optimal distribution is, in first order approximation, the empirical one smoothed over the set of shapes with a Gaussian kernel N ε . Moreover, up to renormalization, g = (Id -ε∆) g 0 coincides with the D T emp of the previous paragraph. The inner products (P i ) which suit g = (g i ) the best in the sense of proposition 1 are precisely the ones obtained in proposition 3 and thus the ones in section 3.1. This is consequently another validation of our approach.

Conclusion

We proposed an approach to learn shape metrics from small training sets of highly-varying shapes, particularly suited to video analysis. The structure we build on sets of shapes relies on deformations and transport, on the contrary to distance-based methods, and allows the consideration of non-dense sample sets. We compute pairwise matchings between close shapes with possibly different topologies, transport deformations with reliability weights, and estimate smooth shape metrics in the whole training set. Thus, we generalize statistical approaches based on deformations, to the case of shape datasets with high variability, where the notion of mean pattern is not relevant anymore.

We studied several ways to estimate metrics, to propose criteria on metrics. We showed that the metric computed in our approach is the optimal one for these criteria, because of a link between Kullback-Leibler divergence and PCA.

We emphasized the new perspectives in segmentation or learning based on shapes, offered by such a transport-based structure. We showed how the metric learned can be turned into a shape matching prior. We also pointed out how to learn all notions (matching, transport) with a second pass, whose completed implementation remains future work.

Figure 1 .

 1 Figure 1. (Left) Example of topological change: hand in the pocket.(Middle) Matching energy explanation. We search for small linking vectors f (s), and for small difference between the two green vectors, or between the red ones, to ensure spatial coherence. These two differences have same norm: ∂f /∂s . (Right) Imposing similar sampling rates limits precision and regularity.

→

Figure 2 .

 2 Figure 2. Colors on one shape, propagated via transports to other shapes from a same video sequence (see text for details).

Figure 4 .

 4 Figure 4. Waving and changing posture: set of 30 hand-segmented shapes from a video. (Left) Transport path between two shapes, (middle) correspondence flow to another shape. A waving sign can be transported to a different posture. (Right) Mean (in black)of the direct matchings from a shape (in red) towards its 5 and 10 nearest neighbors, respectively. Despite small neighborhood size, the mean is irrelevant. Neighborhoods made of reliable matchings are small and cannot include gestures observed at other postures.


Guillaume Charpiat
email: gcharpia@sophia.inria.fr

Matthias Hofmann

Bernhard Schölkopf

Automatic Image Colorization via Multimodal Predictions

We aim to color greyscale images automatically, without any manual intervention. The color proposition could then be interactively corrected by user-provided color landmarks if necessary. Automatic colorization is nontrivial since there is usually no one-to-one correspondence between color and local texture. The contribution of our framework is that we deal directly with multimodality and estimate, for each pixel of the image to be colored, the probability distribution of all possible colors, instead of choosing the most probable color at the local level. We also predict the expected variation of color at each pixel, thus defining a nonuniform spatial coherency criterion. We then use graph cuts to maximize the probability of the whole colored image at the global level. We work in the L-a-b color space in order to approximate the human perception of distances between colors, and we use machine learning tools to extract as much information as possible from a dataset of colored examples. The resulting algorithm is fast, designed to be more robust to texture noise, and is above all able to deal with ambiguity, in contrary to previous approaches.

Introduction

Automatic image colorization consists in adding colors to a new greyscale image without any user intervention. The problem, stated like this, is ill-posed, in the sense that one cannot guess the colors to assign to a greyscale image without any prior knowledge. Indeed, many objects can have different colors: not only artificial, plastic objects can have random colors, but natural objects like tree leaves can have various nuances of green and turn brown during autumn, without a significant change of shape.

The color prior most often considered in the literature is the user: many approaches consist in letting the user determine the color of some areas and then in extending this information to the whole image, either by pre-computing a segmentation of the image into (hopefully) homogeneous color regions, or by spreading color flows from the user-defined color points. This last method supposes a definition of the difficulty of the color flow to go through each pixel, usually estimated as a simple function of local greyscale intensity variations, as in Levin et al. [START_REF] Levin | Colorization using optimization[END_REF], in Yatziv and Sapiro [START_REF] Yatziv | Fast image and video colorization using chrominance blending[END_REF], or in Horiuchi [START_REF] Horiuchi | Colorization algorithm using probabilistic relaxation[END_REF], or by predefined thresholds to detect color edges [START_REF] Takahama | Improvement on colorization accuracy by partitioning algorithm in cielab color space[END_REF]. However, we have noticed, for instance, that the simple, efficient framework by Levin et al. cannot deal with texture examples such as in figure 1, whereas simple oriented texture features such as Gabor filters would have obviously solved the problem 4 . This is the reason why we need to consider texture descriptors. More generally, each hand-made criterion for edge estimation has its drawbacks, and therefore we will learn the probability of color change at each pixel instead of setting a hand-made criterion. User-based approaches present the advantage that the user can interact, add more color points if needed, until a satisfying result is reached, or even place color points strategically in order to give indirect information on the location of color boundaries. Our method can easily be adapted to incorporate such user-provided color information. The task of providing first a fully automatic colorization of the image, before a possible user intervention if necessary, is however much harder. Some recent attempts of predicting the colors gave mixed results. For instance the problem has been studied by Welsh et al. [START_REF] Welsh | Transferring color to greyscale images[END_REF]; it is also one of the applications presented by Hertzmann et al. in [START_REF] Hertzmann | Image analogies[END_REF]. However the framework developed in both cases is not mathematically expressed, in particular it is not clear whether an energy is minimized, and the results shown seem to deal with only a few colors, with many small artifacts probably due to the lack of a suitable spatial coherency criterion. An other notable study is by Irony et al. [START_REF] Irony | Colorization by Example[END_REF]: it consists in finding a few points in the image where a color prediction algorithm reaches the highest confidence, and then in applying Levin's approach as if these points were given by the user. Their approach of color prediction is based on a learning set of colored images, partially segmented by the user into regions. The new image to be colored is then automatically segmented into locally homogeneous regions whose texture is similar to one of the colored regions previously observed, and the colors are transferred. Their method reduces the effort required of the user but still requires a manual pre-processing step. To avoid this, they proposed to segment the training images automatically into regions of homogeneous texture, but fully automatic segmentation (based on texture or not) is known to be a very hard problem. In our approach we will not rely on an automatic segmentation of the training or test images but we will build on a more robust ground.

Irony's method also brings more spatial coherency than previous approaches, but the way coherency is achieved is still very local since it can be seen as a onepixel-radius filter, and furthermore it relies partly on an automatic segmentation of the new image.

But above all, the latter method, as well as all other former methods in the history of colorization (to the knowledge of the authors), cannot deal with the case where an ambiguity concerning the colors to assign can be resolved only at the global level. The point is that local predictions based on texture are most often very noisy and not reliable, so that information needs to be integrated over large regions to become significant. Similarly, the performance of an algorithm based on texture classification, such as Irony's one, would drop dramatically with the number of possible texture classes, so that there is a real need for robustness against texture misclassification or noise. In contrast to previous approaches, we will avoid to rely on very-local texture-based classification or segmentation and we will focus on more global approaches.

The color assignment ambiguity also happens when the shape of objects is relevant to determine the color of the whole object. More generally, it appears that boundaries of objects contain much not-immediately-usable information, such as the presence of edges in the color space, and also contain significant details which can help the identification of the whole object, so that the colorization problem cannot be solved at the local level of pixels. In this article we try to make use of all the information available, without neglecting any low probability at the local level which could make sense at the global level.

Another source for prior information is motion and time coherency in the case of video sequences to be colored [START_REF] Levin | Colorization using optimization[END_REF]. Though our framework can easily be extended to the film case and also benefit from this information, we will deal only with still images in this article.

First we present, in section 2, the model we chose for the color space, as well as the model for local greyscale texture. Then, in section 3, we state the problem of automatic image colorization in machine learning terms, explain why the naive approach, which consists in trying to predict directly color from texture, performs poorly, and we show how to solve the problem by learning multimodal probability distributions of colors, which allows the consideration of different possible colorizations at the local level. We also show how to take into account user-provided color landmarks when necessary. We start section 4 with a way to learn how likely color variations are at each pixel of the image to color. This defines a spatial coherency criterion which we use to express the whole problem mathematically. We solve a discrete version of it via graph cuts, whose result is already very good, and refine it to solve the original continuous problem.

Model for colors and greyscale texture

We first model the basic quantities to be considered in the image colorization problem.

Let I denote a greyscale image to be colored, p the location of one particular pixel, and C a proposition of colorization, that is to say an image of same size as I but whose pixel values C(p) are in the standard RGB color space. Since the greyscale information is already given by I(p), we should add restrictions on C(p) so that computing the greyscale intensity of C(p) should give I(p) back. Thus the dimension of the color space to be explored is intrinsically 2 rather than 3.

We present in this section the model chosen for the color space, our way to discretize it for further purposes, and how to express continuous probability distributions of colors out of such a discretization. We also present the feature space used for the description of greyscale patches.

L-a-b color space

In order to quantify how similar or how different two colors are, we need a metric in the space of colors. Such a metric is also required to associate to any color a corresponding grey level, i.e. the closest unsaturated color. This is also at the core of the color coherency problem: an object with a uniform reflectance will show different colors in its illuminated and shadowed parts since they have different grey levels, so that we need a way to define robustly colors against changes of lightness, that is to consider how colors are expected to vary as a function of the grey level, i.e. how to project a dark color onto the subset of all colors who share a particular brighter grey level.

The psychophysical L-a-b color space was historically designed so that the Euclidean distance between the coordinates of any colors in this space approximates as well as possible the human perception of distances between colors. The transformation from standard RGB colors to L-a-b consists in first applying gamma correction, followed by a linear function in order to obtain the XY Z color space, and then by another highly non-linear application which is basically a linear combination of the cubic roots of the coordinates in XY Z. We refer to Lindbloom's website 5 for more details on color spaces or for the exact formulas. The L-a-b space has three coordinates: L expresses the luminance, or lightness, and is consequently the greyscale axis, whereas a and b stand for two orthogonal color axes.

We choose the L-a-b space to represent colors since its underlying metric has been designed to express color coherency. In the following, by grey level and 2D color we will mean respectively L and (a, b).

With the previous notations, since we already know the grey level I(p) of the color C(p) to predict at pixel p, we search only for the remaining 2D color, denoted by ab(p).

Discretization of the color space

This subsection contains more technical details and can be skipped in the first reading. In section 3 we will temporarily need a discretization of the 2D color space. Instead of setting a regular grid, we define a discretization adapted to the color images given as examples so that each color bin will contain approximately the same number of observed pixels with this color. Indeed, some entire zones of the color space are useless and we prefer to allocate more color bins where the density of observed colors is higher, so that we can have more nuances where it makes statistical sense. Figure 2 shows the densities of colors corresponding to some images, as well as the discretizations in 73 bins resulting from these densities. We colored each color bin by the average color of the points in the bin. To obtain these discretizations, we used a polar coordinate system in ab and cut recursively color bins with highest numbers of points at their average color into 4 parts. Discussing the precise discretization algorithm is not relevant here provided it makes sense statistically; we could have used K-means for instance.

Further, in section 4, we will need to express densities of points over the whole plane ab, based on the densities computed in each color bin. In order to interpolate continuously the information given by each color bin i, we place Gaussian functions on the average color µ i of each bin, with a standard deviation proportional to the statistically observed standard deviation σ i of points of each color bin (see last column of figure 2). Then we interpolate the original density d(i) to any point x in the ab plane by:

d G (x) = i 1 π(ασ i ) 2 e - (x-µ i ) 2 2(ασ i ) 2 d(i)
We found experimentally that considering a factor α ≈ 2 improved significantly the distribution aspect. Cross-validation could be used to determine the optimal α for a given learning set.

Greyscale patches and features

The grey level of one pixel is clearly not informative enough to decide which 2D color we should assign to it. The hope is that texture and local context will give more clues for such a task (see section 3). In order to extract as much information as possible to describe local neighborhoods of pixels in the greyscale image, we compute SURF descriptors [START_REF] Bay | Surf: Speeded up robust features[END_REF] at three different scales for each pixel. This leads to a vector of 192 features per pixel. In order to reduce the number of features and to condense the relevant information, we apply Principal Component Analysis (PCA) and keep the first 27 eigenvectors, to which we add, as supplemental components, the pixel grey level as well as two biologically inspired features: a weighted standard deviation of the intensity in a 5 × 5 neighborhood (whose meaning is close to the norm of the gradient), and a smooth version of its Laplacian. We will refer to this 30-dimensioned vector, computed at each pixel, as local description in the following, and denote it by v or w.

Color prediction

Now that we have modelled colors as well as local descriptions of greyscale images, we can start stating the image colorization problem. Given a set of examples of color images, and a new greyscale image I to be colored, we would like to extract knowledge from the learning set to predict colors C(p) for the new image.

Need for multimodality

One could state this problem in simple machine learning terms: learn the function which associates to any local description of greyscale patches the right color to assign to the center pixel of the patch. This could be achieved by kernel regression tools such as Support Vector Regression (SVR) or Gaussian Processes [START_REF] Schölkopf | Learning with Kernels: Support Vector Machines, Regularization, Optimization, and Beyond[END_REF].

There is an intuitive reason why this would perform poorly. Many objects can have different colors, for instance balloons at a fair could be green, red, blue, etc., so that even if the task of recognizing a balloon was easy and that we knew that we should use colors from balloons to color the new one, a regression would recommend the use of the average value of the observed balloons, i.e. grey. The problem is however not specific to objects of the same class. Local descriptions of greyscale patches of skin or sky are very similar, so that learning from images including both would recommend to color skin and sky with purple, without considering the fact that this average value is never probable.

We therefore need a way to deal with multi-modality, i.e. to predict different colors if needed, or more exactly, to predict at each pixel the probability of every possible color. This is in fact the conditional probability of colors knowing the local description of the greyscale patch around the pixel considered. We will also be interested in the confidence in these predictions in order to know whether some predictions are more or less reliable than others.

Probability distributions as density estimations

The conditional probability of the color c i at pixel p knowing the local description v of its greyscale neighborhood can be expressed as the fraction, amongst colored examples e j = (w j , c(j)) whose local description w j is similar to v, of those whose observed color c(j) is in the same color bin B i . We thus have to estimate densities of points in the feature space of grey patches. This can be accomplished with a Gaussian Parzen window estimator:

p(c i |v) = {j : c(j)∈Bi} k(w j , v) j k(w j , v)
where k(w j , v) = e -(wj -v) 2 /2σ 2 is the Gaussian kernel. The best value for the standard deviation σ can be estimated by cross-validation on the densities. With this framework we can express how reliable the probability estimation is: its confidence depends directly on the density of examples around v, since an estimation far from the clouds of observed points loses signification. Thus, the confidence in a probability prediction is the density in the feature space itself:

p(v) ∝ j k(w j , v)
In practice, for each pixel p, we compute the local description v(p), but we do not need to compute the similarities k(v, w j ) to all examples in the learning set: in order to save computational time, we only search for the K-nearest neighbors of v in the learning set, with K sufficiently large (as a function of the σ chosen), and estimate the Parzen densities based on these K points. In practice we choose K = 500, and thanks to fast nearest neighbor search techniques such as kD-tree6 , the time needed to compute all predictions for all pixels of a 50 × 50 image is only 10 seconds (for a learning set of hundreds of thousands of patches) and this scales linearly with the number of test pixels. Note that we could also have used sparse kernel techniques such as SVR to estimate, for each color bin, a regression between the local descriptions and the probability of falling into the color bin. We refer more generally to [START_REF] Schölkopf | Learning with Kernels: Support Vector Machines, Regularization, Optimization, and Beyond[END_REF] for details and discussions about kernel methods.

User-provided color landmarks

We can easily consider user-provided information such as the color c at pixel p in order to modify a colorization obtained automatically. We set p(c|p) = 1 and set the confidence p(p) to a very large value. Consequently our optimization framework is still usable for further interactive colorization. A re-colorization with new user-provided color landmarks does not require the re-estimation of color probabilities, and therefore lasts only a fraction of second (see next part).

Spatial coherency

For each pixel of a new greyscale image, we are now able to estimate the probability distribution of all possible colors (within a big finite set of colors since we discretized the color space into bins). The interest in such computation is that, if we add a spatial coherency criterion, a pixel will be influenced by its neighbors, and the choice of the best color to assign will be done accordingly to the probability distributions in the neighborhood. Since all pixels are linked by neighborhoods, even if not directly, they all interact with each other, so that the solution has to be computed globally. Indeed it may happen that, in some regions that are supposed to be homogeneous, a few different colors may seem to be the most probable ones at a local level, but that the winning color at the scale of the region is different, because in spite of its only second rank probability at the local level, it ensures a good probability everywhere in the whole region. The opposite may also happen: to flip a whole such region to a color, it may be sufficient that this color is considered as extremely probable at a few points with high confidence. The problem is consequently not trivial, and the issue is to find a global solution. In this section we first learn a spatial coherency criterion, then find a good solution to the whole problem with the help of graph cuts.

Local color variation prediction

Instead of picking randomly a prior for spatial coherence, based either on detection of edges, or on the Laplacian of the intensity, or on a pre-estimated complete segmentation, we learn directly how likely it is to observe a color variation at a pixel knowing the local description of its greyscale neighborhood, based on a learning set of real color images. The technique is similar to the one detailed in the previous section. For each example w j of colored patch, we compute the norm g j of the gradient of the 2D color (in the L-a-b space) at the center of the patch. The expected color variation g(v) at the center of a new greyscale patch v is then:

g(v) = j k(w j , v) g j j k(w j , v)
.

Thus we now have priors both on the colors and on the color variations.

Global coherency via graph cuts

The graph cut, or max flow, algorithm is a minimization technique widely used in computer vision [START_REF] Boykov | An experimental comparison of min-cut/max-flow algorithms for energy minimization in vision[END_REF][START_REF] Kolmogorov | What energy functions can be minimized via graph cuts?[END_REF] because of its suitability for many image processing problems, because of its guarantee to find a good local minimum, and because of its speed. In the multi-label case with α-expansion [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF], it can be applied to all energies of the form i V i (x i ) + i∼j D i,j (x i , x j ) where x i are the unknown variables, with possible values in a finite set L of labels, where the V i are any functions, and where D i,j are any pair-wise interaction terms with the restriction that each D i,j (•, •) should be a metric on L.

The reason why we temporarily discretized the color space in section 2 was to be able to use this technique. We formulate the image colorization problem as an optimization problem on the following energy:

p V p (c(p)) + λ p∼q |c(p) -c(q)| Lab g p,q (1)
where

V p (c(p)) = -log p v(p) p c(p)|v(p)
is the cost of choosing color c(p) for pixel p (whose neighboring texture is described by v(p))) and where g p,q = 2 g(v(p)) -1 + g(v(q)) -1 -1 is the harmonic mean of the estimated color variation at pixels p and q. An 8-neighborhood is considered for the interaction term, and p ∼ q means that p and q are neighbors.

The term V p penalizes colors which are not probable at the local level according to the probability distributions obtained in section 3, with the strength depending on the confidence in the predictions. The interaction term between pixels penalizes color variation where it is not expected, according to the variations predicted in the previous paragraph.

We use the graph cut package7 provided by [START_REF] Szeliski | A comparative study of energy minimization methods for markov random fields[END_REF]. The solution for a 50 × 50 image and 73 possible colors is obtained by graph cuts in a fraction of second and is generally already very satisfying. The computation time scales approximately quadratically with the size of the image, which is still fast, and the algorithm performs well even on massively down-scale versions of the image, so that a good initial clue can still be given quickly for very big images too. The computational costs compete with those of the fastest colorization techniques [START_REF] Blasi | Fast Colorization of Gray Images[END_REF] while bringing much more spatial coherency.

Refinement in the continuous color space

We can now go back to the initial problem in the continuous space of colors. We interpolate probability distributions p(c i |v(p)) estimated at each pixel p for each color bin i, to the whole space of colors with the technique described in section 2, so that V p (c) is now defined for any color c and not only for color bins. The energy (1) can consequently be minimized in the continuous space of colors. We start from the solution obtained by graph cuts and refine it with a gradient descent. This refinement step will generally not introduce huge changes like flipping of whole regions, but will bring more nuances. Fig. 3. Da Vinci case: Mona Lisa colored with Madonna. The border is not colored because of the window size needed for SURF descriptors. Second line: color variation predicted (white stands for homogeneity and black for color edge); most probable color at the local level; 2D color chosen by graph cuts. Note that previous algorithms could not deal with regions such as the neck or the forehead, where blue is the most probable color at the local level because greyscale skin looks like sky. Surroundings of these regions and lower-probability colors are decisive for the final color choice.

Experiments

We now show results of automatic colorization. In figure 3 we colored a famous painting by Leonardo Da Vinci with another painting of his. The paintings are significantly different and textures are relatively dissimilar. The prediction of color variation performs well and helps much to determine the boundaries of homogeneous color regions. The multimodality framework proves extremely useful in areas such as Mona Lisa's forehead or neck where the texture of skin can be easily mistaken with the texture of sky at the local level. Without our global optimization framework, several entire skin regions would be colored in blue, disregarding the fact that skin color is a second probable possibility of colorization for these areas, which makes sense at the global level since they are surrounded by skin-colored areas, with low probability of edges. We insist on the fact that the input of previous texture-based approaches is very similar to the "most probable color" prediction (second line, middle image), whereas we consider the probabilities of all possible colors at all pixels. This means that, given a certain quality of texture descriptors, we handle much more information.

In figure 4 we perform similar experiments with photographs of landscapes. The effect of the refinement step can be observed in the sky, where nuances of blue vary more smoothly. 3, plus last image: colors obtained after refinement step. Note that the sky is more homogeneous, the color gradients in the sky are smoother than when obtained directly by graph cuts (previous image).

We compare our method with the one by Irony et al., on their own example [START_REF] Irony | Colorization by Example[END_REF] in figure 5; the task is easier and results are similar. The boundaries of our color regions even fit better to the zebra contour. Grass areas near the zebra are colored according to the grass observed at similar locations around the zebra in the learning image. This bias should disappear with bigger training sets since the color of the background would become independent of zebra's presence. We would have liked to compare results on a wider scale, in particular on our difficult examples, but we face the problems of source code availability and objective comparison criteria.

In figure 6 we consider a very difficult task: the one of coloring an image from a Charlie Chaplin movie, with many different objects and textures, such as a brick wall, a door, a dog, a head, hands, a loose suit... Because of the number of objects, and because of their particular arrangement, it is unlikely to find a single color image with a similar scene that we would use as a learning image. Thus we consider a small set of three different images, each of which shares a partial similarity with Charlie Chaplin's image. The underlying difficulty is that each learning image also contains parts which should not be re-used in this target image. The result is however promising, considering the learning set. Dealing with bigger learning datasets should slow the process only logarithmically, during the kD-tree search.

Discussion

Our contribution is both theoretical and practical.

We proposed a new approach for automatic image colorization, which does not require any intervention by the user, except from the choice of relatively similar color images. We stated the problem mathematically as an optimization problem with an explicit energy to minimize. Since it deals with multimodal probability distributions until the final step, this approach makes better use of the information that can be extracted from images by machine learning tools. The fact that we solve directly the problem at the global level, with the help of graph cuts, makes the framework more robust to noise and local prediction errors. It also makes it possible to resolve large scale ambiguities which previous approaches could not. This multi-modality framework is not specific to image colorization and could be re-used in any prediction task on images.

We also proposed a way to learn and predict where color variations are probable and how important they are, instead of choosing a spatial coherency criterion Fig. 6. Charlie Chaplin frame colored by three different images. Second line: prediction of color variations, most probable color at the local level, and final result. This example is particularly difficult because many different textures and colors are involved and because there does not exist color images very similar to Charlie Chaplin's one. We have chosen three different images, each of which shares a partial similarity with the target image. In spite of the complexity, the prediction of color edges or of homogeneous regions remains significant. The brick wall, the door, the head and the hands are globally well colored. The large trousers are not in the learning set; the mistakes in the colors of Charlie Chaplin's dog are probably due to the blue reflections on the dog in the learning image and to the light brown of its head. Both would probably benefit from larger training sets.

by hand, and this performs quite well (see figures 3 and 6 again). It relies on many features and adapts itself to the current training set.

Our colorization approach outperforms [START_REF] Welsh | Transferring color to greyscale images[END_REF] and [START_REF] Hertzmann | Image analogies[END_REF], whose examples contain only few colors and lack spatial coherency. Our process requires less or similar intervention than [START_REF] Irony | Colorization by Example[END_REF] but can handle more ambiguous cases and more texture noise. The computational time needed is also very low, enabling real-time user interaction since a re-colorization with new color landmarks lasts a fraction of second.

The automatic colorization results should not be compared to those obtained by user-helped approaches since we do not benefit from such decisive information. Nevertheless, even with such an handicap, our results clearly compete with the state of the art of colorization with few user scribbles.

Fig. 1 .

 1 Fig. 1. Failure of standard colorization algorithms in presence of texture. Left: manual initialization; right: result of Levin et al.'s approach. In spite of the general efficiency of their simple method (based on the mean and standard deviation of local intensity neighborhoods), texture remains difficult to deal with. Hence the need of texture descriptors and of learning edges from color examples.

Fig. 2 .

 2 Fig. 2. Examples of color spectra and associated discretizations. For each line, from left to right: color image; corresponding 2D colors; location of the observed 2D colors in the ab-plane (a red dot for each pixel) and the computed discretization in color bins; color bins filled with their average color; continuous extrapolation: influence zones of each color bin in the ab-plane (each bin is replaced by a Gaussian, whose center is a black dot; red circles indicate the standard deviation of colors within the color bin, blue ones are three times bigger).

Fig. 4 .

 4 Fig.4. Landscape example. Same display as in figure3, plus last image: colors obtained after refinement step. Note that the sky is more homogeneous, the color gradients in the sky are smoother than when obtained directly by graph cuts (previous image).

Fig. 5 .

 5 Fig. 5. Comparable results with Irony et al. on their own example[START_REF] Irony | Colorization by Example[END_REF]. First line: our result. Second line, left: our predicted 2D colors; right: Irony et al.'s result with the assumption that it is a binary classification problem. Note that this example they chose could be solved easily without texture consideration since there is a simple correspondence between grey level intensity and color. Concerning our result, the color variations in the grass around the zebra are probably due to the influence of the grass color for similar patches in the learning image, and this problem should disappear with a larger training set. Note however, if you zoom in, that the contour of the zebra matches exactly boundaries of 2D color regions (purple and pink), whereas in Irony et al.'s result, a several pixel wide band of grass around the zebra's legs and abdomen are colored as if they were part of the zebra. The moral is that we should not try to learn colors from a single overfitted, carefully-chosen training image but rather check how the colorization scales with the number of training images and with their non-perfect suitability.

Their code is available at http://www.cs.huji.ac.il/~weiss/Colorization/.

http://brucelindbloom.com/

We use, without particular optimization, the TSTOOL package available at http://www.physik3.gwdg.de/tstool/.

available at http://vision.middlebury.edu/MRF/code/

Acknowledgments

We would like to thank Jason Farquhar, Peter Gehler, Matthew Blaschko and Christoph Lampert for very fruitful discussions.