À Aurélie

Esteban Anaé

Keywords: Generalized Wasserstein distance, image processing, image interpolation, color transfer, segmentation, (non)convex optimization Distance de Wasserstein généralisée, interpolation d'images, transfert de couleur, segmentation, optimisation

Optimal Transport is a well developed mathematical theory that defines robust metrics between probability distributions. The computation of optimal displacements between densities through the associated transport map makes this theory mainstream in several applicative fields. For image processing applications, the transport map can for instance be used to compute geodesics between images or to transfer characteristics of one image to another. In this context, it is of major interest to preserve the nature of the observed objects, in order to synthesize images that are physically and visually plausible. In this document, generalized Optimal Transport distances including relaxation and regularization are considered to improve the modeling of image processing problems. New models and algorithms within the continuous and discrete formulations of Optimal Transport are presented. With the continuous setting, the integration of physical regularization of the transport plan makes possible the interpolation of ocean images containing complex structures. In the discrete setting, the regularization of the transport plan is considered for color transfer between images. Convex and non-convex models are proposed to define automatic methods that adapt the proportion of colors required to synthesize visually plausible images. These methods are extended to the computation of barycenters to deal with the color normalization of multiple images. Finally, the entropy regularization of discrete Optimal Transport is used for image segmentation. A fast and convex model is designed to segment images, while respecting global color distribution constraints.

Conclusion of Part I

Chapter 1 Introduction

This habilitation manuscript presents a selected subset of the research activities realized in the period 2011-2015. It focuses on the use of Optimal Transport for image processing purposes. These works have been initiated through the supervision of two PhD thesis and one post-doctorate.

In this introduction, the applicative and scientific motivations of the realized works are presented. The main limitations of the use of Optimal Transport for Image Processing are then exposed. The related contributions are finally briefly presented, before detailing the organization of the document.

Motivations

Computer Vision Image and video contents are omnipresent in the world today, both in our personal and professional lives. Technology is in constant evolution and the amount of data is everyday dramatically increasing. Market data indicates that the number of acquisition mobile devices is going to considerably increase. In particular, a spectacular growth in video users for the years to come is expected. In 2016, more than 1.6 billion electronic devices capable of recording and sharing pictures and videos will be used all around the world. Billions of raw images and videos are consequently diffused on the Internet. For instance, approximately 2 millions of pictures are uploaded on Flickr on a daily basis and the volume of videos uploaded to Youtube is gigantic, and constantly increasing (300 hours per minute). However, most of this video content is not exploitable because it has not been properly edited.

Due to the diversity of applications in which images are involved, the type of data is completely heterogeneous: color images, 3D images or even multidimensional images (including transparency, depth, etc), animated sequence and videos, multispectral images, depth images, light-field, etc. Additionally, pictures are often captured at different resolutions and under different lighting conditions.

The need for generic methods to deal with this huge amount of data is increasing. The development of methods for a given application implies the control of the statistics of the images targeted by this application. An example is presented in Figure 1.1. Raw images diffused on the Internet are used for 3D reconstruction of buildings (using e.g. geolocalized Flickr images) or virtual synthesis of 2D views of cities (Google street-view, Institut National de l'Information Géographique et Forestière, IGN). For these applications, image harmonization is a key step to merge heterogeneous data with different color statistics. In order to circumvent this challenging task, the current methods only select a small subset of available images that are sufficiently close in terms of color and illuminations. As a consequence, most of existing data is not even used, whereas it could lead to significant improvements. In order to use more images in this context or simply when one looks for pictures on the Internet to illustrate professional or personal presentations, the huge quantity of images needs to be efficiently indexed and retrieved. This requires very fast and efficient comparison tools between descriptors of image statistics, which are also of interest for object segmentation, editing or harmonization.

Histograms are popularly used in image processing, computer vision and machine learning to represent complex visual objects. To that end, relevant features such as color, contours, orientations or textures are first extracted from images. Histograms are built from the quantization of the space of features into discrete bins. One can further define normalized histograms that describe the frequencies of each of the observed features. The problems of retrieving a similar image in a database or segmenting a particular object can thus be recasted as a problem of histograms comparison. In this context, even if adapted metrics can be learned by neural networks for specific application, generic robust comparison tools are still needed for a broad use.

While a lot of research effort is dedicated to the enhancement of image descriptors (by adding additional information such as text or context) in a "big data" perspective, the common metrics used for statistics comparison are still not developed enough to deal with various perturbation effects (such as histogram quantization, shifts, deformation, etc) or large amount of data outliers. Namely, spatial information is not efficiently considered as it is mainly used to enhance the feature description [START_REF] Sharma | Expanded parts model for human attribute and action recognition in still images[END_REF]. Including spatial information into the distance defined on the feature space is an interesting alternative to improve the quality of the results.

The Optimal Transport framework is the proper way to compare statistic distributions. In contrast to most distances from information theory, it takes into account the spatial location of the density modes [START_REF] Villani | Topics in Optimal Transportation[END_REF]. By incorporating in its definition a ground metric between the features themselves, the Optimal Transport distance (that corresponds to the Wasserstein distance when the ground metric is the quadratic Euclidean distance) can compare sparse histograms even if their supports do not overlap signif-icantly. Additionally, Optimal Transport provides a warping (the so-called transport plan) between histograms which can be used to perform image editing such as color transfer.

Optimal Transport framework has been shown to produce robust state of the art results for the comparison of statistical descriptors through the so-called Earth Mover's Distance [START_REF] Rubner | A metric for distributions with applications to image databases[END_REF]. Various works on color image retrieval [START_REF] Pele | A linear time histogram metric for improved sift matching[END_REF] or artistic image indexation [START_REF] Hurtut | Adaptive image retrieval based on the spatial organization of colors[END_REF] have shown on small datasets (i.e. thousands of images) that Optimal Transport framework is intrinsically designed to carry this spatial information. However, unless considering some crude approximations and/or simple ground metrics, it involves prohibitive computational costs. Though Optimal Transport would lead to a significant gap in performance improvement, this computational limitation is clearly an obstacle to its study and large scale deployment, when considering the dimension of image databases or the data itself (multi-spectral images, patch representation, tensor-field imaging, etc).

Geosciences Since the late seventies, many satellites have been launched to improve our knowledge of the atmosphere and the oceans. Geostationary satellites provide, among other data, photographic images of the earth system. Sequences of such images show the dynamical evolution of identified meteorological or oceanic "objects": fronts, clouds, eddies, vortices, etc. The human vision can easily detect the dynamics in this kind of image sequences and it clearly has a strong predictive potential. This aspect is favoured by the fact that image data, contrary to many other measurements, is dense in space and time. Indeed the spatial resolution of current METEOSAT satellites is close to one kilometer and they produce an image every 15 minutes. This frequency will be improved up to one every 10 minutes (and even every 2.5 minutes for Europe only) with the upcoming third satellite generation. It implies a huge quantity of information which can be seen as an asset but also induces difficulties for the assimilation system that has to cope with such amount of data.

Satellite data is currently used in operational systems for calibrating Numerical Weather Prediction models, mainly through the assimilation of the radiance measured by the satellite at each pixel of the image. As illustrated in Figure 1.2, satellite images are related to physical quantities such as surface temperature, sea surface height, cloud pressure, chlorophyll concentration, etc. In practice only a tiny percentage (about 3 -5%) of total satellite (from polar orbiting and geostationary) data is used in operational systems, while being given low confidence with respect to synoptic data, collected by stations, aircrafts, radiosounding, ballons or drifters. Considering the cost of satellite observing systems (the cost of the launch of the Meteosat Third Generation is estimated to about 2.5 billion Euros) and of the infrastructures required for the collection of the data itself, improving their impact on forecasting systems is an important topic of research.

As a consequence, there is a need to provide pertinent distances to interpolate between the observed satellite images and the images of physical variables provided by the numerical models. The spatial information contained in satellite images is currently not taken into account, as the classical "pixel-to pixel" L 2 distance used in data assimilation is not adapted to deal with bad localization of structures such as clouds, vortex or fronts [START_REF] Titaud | Assimilation of Image Sequences in Numerical Models[END_REF], [START_REF] Chabot | Accounting for observation errors in image data assimilation[END_REF]. If the same structure is contained in both images at different locations, an interpolation should deal with the spatial shift. While a L 2 interpolation would create two structures, Optimal Transport is intrinsically designed to correctly interpolate the position of the structure along the geodesic in the Wasserstein space.

Considering images as densities, interpolation between images can be realized with the Optimal Transport map [START_REF] Haker | Optimal mass transport for registration and warping[END_REF]. It is thus an interesting framework in order to propose new metrics for data assimilation [START_REF] Papadakis | Assimilation d'images et de structures[END_REF]. In this context, an important challenge is to model the topology of the considered domain (coasts or islands in oceanography) for estimating image interpolations that have a real physical meaning. In this context, even if Optical Flow methods can be considered for data assimilation to compute discrepancies between atmospheric images [START_REF] Corpetti | Pressure image assimilation for atmospheric motion estimation[END_REF], they can not deal with complex domains in oceanography [START_REF] Beyou | Weighted ensemble transform kalman filter for image assimilation[END_REF].

Optimal Transport Optimal Transport is a well developed mathematical theory that defines a family of metrics between probability distributions [START_REF] Gangbo | The geometry of optimal transportation[END_REF][START_REF] Villani | Topics in Optimal Transportation[END_REF]. These metrics measure the amplitude of an optimal displacement according to a so-called ground cost defined on the space supporting the distributions. The resulting distance is referred to as the Wasserstein distance in the case of L p ground costs C(x, y) = ||x -y|| p . As shown in Figure 1.3, it measures the minimal effort required for filling a "remblai" (-µ 1 ) with a "déblai" (µ 0 ), i.e. transporting one distribution to another. The geometric nature of Optimal Transport, as well as the ability to compute optimal displacements between densities through the corresponding transport map T , make this theory progressively mainstream in several applicative fields. Figure 1.3 -Illustration of Optimal Transport between distributions µ 0 and µ 1 as introduced in [START_REF] Monge | Mémoire sur la théorie des déblais et de remblais[END_REF] with the "déblais" and "remblais" problem . Interpolations (µ t ) t∈[0;1] between µ 0 and µ 1 can be compted with the transport map T .

Early successes of applications of Optimal Transport were mostly theoretical, such as the study of shape recognition [START_REF] Gangbo | Shape recognition via wasserstein distance[END_REF], the derivation of functional inequalities [START_REF] Cordero-Erausquin | A mass-transportation approach to sharp Sobolev and Gagliardo-Nirenberg inequalities[END_REF], the construction of solutions of non-linear partial differential equations [START_REF] Jordan | The variational formulation of the Fokker-Planck equation[END_REF] or the study of gradient flows in Wasserstein spaces [START_REF] Ambrosio | Gradient flows: in metric spaces and in the space of probability measures[END_REF].

In computer vision, the Wasserstein distance has been shown to outperform other metrics between distributions for machine learning tasks [START_REF] Rubner | A metric for distributions with applications to image databases[END_REF][START_REF] Pele | Fast and robust earth mover's distances[END_REF][START_REF] Cuturi | Fast computation of wasserstein barycenters[END_REF] or image segmentation [START_REF] Ni | Local histogram based segmentation using the wasserstein distance[END_REF][START_REF] Peyré | Wasserstein active contours[END_REF][START_REF] Swoboda | Variational image segmentation and cosegmentation with the wasserstein distance[END_REF].

In image processing, the warping provided by the Optimal Transport has been used for video restoration [START_REF] Delon | Movie and video scale-time equalization application to flicker reduction[END_REF], color transfer [START_REF] Pitié | The linear monge-kantorovitch linear colour mapping for example-based colour transfer[END_REF], texture synthesis [START_REF] Ferradans | Static and dynamic texture mixing using optimal transport[END_REF], optical nanoscopy [START_REF] Brune | 4D imaging in tomography and optical nanoscopy[END_REF] and medical imaging registration [START_REF] Haker | Optimal mass transport for registration and warping[END_REF]. It has also been applied to interpolation in computer graphics [START_REF] Bonneel | Displacement interpolation using Lagrangian mass transport[END_REF][START_REF] Solomon | Convolutional wasserstein distances: Efficient optimal transportation on geometric domains[END_REF] and surface reconstruction in computational geometry [START_REF] Digne | Featurepreserving surface reconstruction and simplification from defect-laden point sets[END_REF]. Optimal Transport is either used to model various physical phenomena, such as for instance in astrophysics [START_REF] Frisch | A reconstruction of the initial conditions of the universe by optimal mass transportation[END_REF] and oceanography [START_REF] Benamou | Weak existence for the semigeostrophic equations formulated as a coupled Monge-Ampere/transport problem[END_REF][START_REF] Feyeux | Application of optimal transport to data assimilation[END_REF].

Limitations of Optimal Transport in Imaging

The numerical resolution of the Optimal Transport problem raises several challenges. This is the main reason why it has been poorly used in image processing until recently. Computing Optimal Transport distance is only an easy task when dealing with one-dimensional densities, such as histograms of grayscale images, or multi-dimensional densities roughly discretized. However, when considering transport between images themselves or when dealing with statistics of color images, the densities are defined on higher dimensional spaces. Hence, the versatility and high quality of Optimal Transport distances come at a price for general image applications: they require the resolution of high dimensional problems that scale with the product of the dimensions of the discretized densities.

In this context, the fluid dynamic formulation of the Optimal Transport problem introduced in [START_REF] Benamou | A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem[END_REF] is an interesting approach for dealing with higher dimensions. The entropic regularization proposed in [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF] has also offered new perspectives as faster algorithms based on Sinkhorn distances can be designed for computing approximate Optimal Transport in larger dimensions. This last approach gives interesting approximations of the transport distance but leads to transport maps that may be too far from the optimal ones for transfer and interpolation purposes.

Hence, generalizations of the Optimal Transport distances are still necessary to adapt to applications involving image processing. For instance, when the transport plan is used to interpolate between densities that have very different "shapes" or contain data outliers, partial models able to relax the constraint of transporting the whole mass are required. In this context, another flaw for the use of Optimal Transport plan is that it is usually highly irregular. For the interpolation between images or histograms of features, dedicated regularizations of the transport plan are needed to preserve either the nature of the objects contained in the images or the modes of the statistics.

Such property is illustrated in Figure 1.4 with an example of color transfer between images. If the exact prescription of color is realized with the Optimal Transport map, visual artifacts may appear. Finally, introducing Optimal Transport distances in more general image problems such as segmentation, involving histogram data fitting and spatial regularization in the image domain, is even more challenging. Spatial regularization can be incorporated through graph or variational modeling. However, graph based approaches are limited to simple bin-to-bin metrics between histograms [START_REF] Rother | Cosegmentation of image pairs by histogram matching-incorporating a global constraint into MRFs[END_REF] and cannot deal with Optimal Transport distances without leading to problems of gigantic size [START_REF] Solomon | Wasserstein propagation for semi-supervised learning[END_REF]. Slow hybrid methods that alternate optimization on graph and Optimal Transport computation should be designed. Variational methods may thus be more adapted to such problems, since the computation of Optimal Transport distances can be integrated into a general model that can directly be optimized.

Contributions and Organization of the document

In this document, some new variational image processing models using Optimal Transport distances are exposed. The main common contribution of the presented works is to propose original modeling of classical problems such as image interpolation, color transfer or image segmentation, while considering recent optimization methods able to solve these problems with "reasonable" computational cost (several seconds to several minutes).

For that purpose, extended Optimal Transport models including relaxation of the mass constraint and regularization of the transport map are introduced. These generalized Optimal Transport models have been designed in both continuous and discrete settings. The continuous formulation is applied to image interpolation, while the discrete one is used for color transfer and image segmentation.

The problems associated to both formulations of Optimal Transport are illustrated in Figure 1.5. In the continuous framework, the images are considered as continuous densities µ 0 and µ 1 defined on a domain Ω. The transport map is in this case a vector field, defined from Ω to Ω, that transfers µ 0 onto µ 1 . The problem can be expressed in a dynamic way through a fluid mechanics formulation. In the discrete setting, normalized histograms of features extracted from images are considered. When computing the static transport between histograms X and Y discretized with M and N bins, an acceptable transport map is a M × N matrix describing a joint distribution P which marginals are the given histograms (i.e. j P i,j = X i and i P i,j = Y j ). Hence P i,j is the mass transferred from X i to Y j . Overview of the manuscript The first part of this document is dedicated to some contributions in the continuous setting. Chapter 2 is based on [START_REF] Papadakis | Optimal transport with proximal splitting[END_REF] and presents an approach for the numerical estimation of the geodesic path between two densities according to the L 2 Wasserstein metric. Generalizations including additional physical constraints proposed during the PhD of R. Hug in [START_REF] Hug | Multiphysics optimal transportation and image interpolation[END_REF] for application to the interpolation of ocean images are next presented in Chapter 3.

The second part focuses on contributions in the discrete setting. Chapter 4 presents the relaxed and regularized formalism proposed in [START_REF] Ferradans | Regularized discrete optimal transport[END_REF] as well as entropic regularization recently introduced in the literature [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF]. Chapters 5 and 6 are respectively dedicated to the application of these tools to color transfer between images and image segmentation. Chapter 5 gathers methods initiated during the Postdoctorate of S. Ferrandans [START_REF] Ferradans | Regularized discrete optimal transport[END_REF][START_REF] Rabin | Non-convex relaxation of optimal transport for color transfer[END_REF], while Chapter 6 reviews some models introduced during the PhD of R. Yıldızoğlu [START_REF] Yıldızoğlu | Active contours without level sets[END_REF][START_REF] Yıldızoğlu | A convex formulation for global histogram based binary segmentation[END_REF][START_REF] Rabin | Convex color image segmentation with optimal transport distances[END_REF].

In order to make the manuscript self-content, standard first-order proximal algorithms dedicated to the optimization of non-smooth functionals used all along the document are finally presented in Appendix A.

Chapter 2 Dynamic Optimal Transport

This chapter is focused on the computation of geodesics for the Optimal Transport metric associated to the L 2 ground cost. It reviews various methods and extends the approach pioneered in [START_REF] Benamou | A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem[END_REF] from the perspective of proximal operator splitting in convex optimization [START_REF] Papadakis | Optimal transport with proximal splitting[END_REF]. It shows the simplicity and efficiency of this method, which can easily be extended beyond the setting of Optimal Transport by considering various convex cost functions that will be presented in Chapter 3.

Estimation of the Transport Map

Let Ω be a convex and bounded domain of R d and ρ 0 , ρ 1 be two non-negative L 1 functions on Ω, of equal mass. We will assume without loss of generality that

Ω ρ 0 (x) dx = Ω ρ 1 (x) dx = 1.
(2.1)

The original formulation of the Optimal Transport problem in [START_REF] Monge | Mémoire sur la théorie des déblais et de remblais[END_REF] corresponds to minimizing the cost for transporting a ρ 0 onto ρ 1 using a map T . In the following, we restrict our exposition to maps

T : Ω ⊂ R d → Ω. A valid transport map T : Ω → Ω then pushes forward the measure ρ 0 (x)dx onto ρ 1 (x) dx, that is T #ρ 0 = ρ 1 or for all bounded set A ⊂ Ω: A ρ 1 (x)dx = T (x)∈A ρ 0 (x)dx. (2.2)
In term of densities, when T is smooth and one-to-one, the constraint

T #ρ 0 = ρ 1 corresponds to ρ 0 (x) = ρ 1 (T (x)) |det(∇T (x))| (2.3)
where ∇T (x) ∈ R d×d is the differential of T at x. This is known as the Jacobian equation. We call T (ρ 0 , ρ 1 ) the set of transport maps that satisfy the constraint (2.3). The L p Kantorovitch-Wassertein distance between ρ 0 and ρ 1 is then defined by

W p (ρ 0 , ρ 1 ) p = inf T ∈T (ρ 0 ,ρ 1 )
|T (x) -x| p ρ 0 (x)dx.

Estimation of the Transport Map

The L p Monge-Kantorovitch problem corresponds to finding a mapping T such that this infimum is achieved. It can be recasted as

W p (ρ 0 , ρ 1 ) p = min T ∈T (ρ 0 ,ρ 1 ) C(x, T (x))ρ 0 (x) dx (2.4)
where C(x 0 , x 1 ) ≥ 0 is the ground cost of transporting x 0 ∈ Ω onto x 1 ∈ Ω.

Let us now detail how estimating the transport map T with Partial Differential Equations.

Optimal Transport and PDEs

The Optimal Transport for the L 2 ground cost has a special structure. It can be shown to be uniquely defined (see e.g. [START_REF] Villani | Topics in Optimal Transportation[END_REF] page 66) and the transport map is the gradient of a convex potential Ψ from Ω to R:

T (x) = ∇Ψ(x).
(2.5)

Hence, one can show that the mass transfer associated to a L 2 ground cost follows straight lines [START_REF] Brenier | Polar factorization and monotone rearrangement of vector-valued functions[END_REF], which can be used for developing specific Lagrangian solvers [START_REF] Iollo | A lagrangian scheme for the solution of the optimal mass transfer problem[END_REF].

Other class of methods relies on the fact that from relations (2.3) and (2.5), the convex function Ψ is solution of the Monge-Ampère equation:

det(D 2 Ψ)ρ 1 (∇Ψ(x)) = ρ 0 (x).
This equation being highly nonlinear, numerical methods to solve the Monge-Kantorovitch problem based on discretization of the Monge-Ampère equation have been investigated [START_REF] Oliker | On the numerical solution of the equation ∂ 2 z ∂x 2 ∂ 2 z ∂y 2 -∂ 2 z ∂x∂y 2 = f and its discretizations[END_REF][START_REF] Loeper | Numerical solution of the Monge-Ampère equation by a Newton's algorithm[END_REF][START_REF] Dean | An augmented Lagrangian approach to the numerical solution of the dirichlet problem for the elliptic Monge-Ampère equation in two dimensions[END_REF][START_REF] Oberman | Wide stencil finite difference schemes for the elliptic Monge-Ampere equation and functions of the eigenvalues of the hessian[END_REF][START_REF] Benamou | Numerical solution of the optimal transportation problem via viscosity solutions for the Monge-Ampère equation[END_REF][START_REF] Benamou | Numerical solution of the optimal transportation problem using the Monge-Ampère equation[END_REF]. A major difficulty in these approaches is to deal with compactly supported densities, which requires a careful handing of the boundary conditions [START_REF] Froese | A numerical method for the elliptic Monge-Ampère equation with transport boundary conditions[END_REF].

Another line of methods iteratively constructs mass preserving mappings converging to the Optimal Transport [START_REF] Angenent | Minimizing flows for the Monge-Kantorovich problem[END_REF][START_REF] Haber | An efficient numerical method for the solution of the l 2 optimal mass transfer problem[END_REF][START_REF] Bouharguane | Numerical solution of the Monge-Kantorovich problem by density lift-up continuation[END_REF]. This explicitly constructs the so-called polar factorization of the initial map, see also [START_REF] Benamou | A domain decomposition method for the polar factorization of vector-valued mappings[END_REF] for a different approach. These PDE's based approaches to the resolution of the Optimal Transport have found several applications, such as image registration [START_REF] Haker | Optimal mass transport for registration and warping[END_REF], density regularization [START_REF] Burger | Regularized regression and density estimation based on optimal transport[END_REF], optical flow [START_REF] Clarysse | Optimal extended optical flow subject to a statistical constraint[END_REF] and grid generation [START_REF] Sulman | Optimal mass transport for higher dimensional adaptive grid generation[END_REF].

A last axis of research consists in using gradient flows where the gradient direction is computed according to the Wasserstein distance. This was initially proposed in [START_REF] Jordan | The variational formulation of the Fokker-Planck equation[END_REF] to build solutions to certain non-linear PDE's. This technique is now being used to design numerical approximation schemes for the solution of these equations, see for instance [START_REF] Carrillo | Numerical simulation of diffusive and aggregation phenomena in nonlinear continuity equations by evolving diffeomorphisms[END_REF][START_REF] During | A gradient flow scheme for nonlinear fourth order equations[END_REF][START_REF] Ferragut | Mixed finite element methods for a class of nonlinear reaction diffusion problems[END_REF][START_REF] Benamou | Discretization of functionals involving the Monge-Amp\'ere operator[END_REF].

All these approaches are nevertheless limited to positive densities, or to non-negative densities with a convex support of the target density ρ 1 , which is a real limitation for image processing applications. The only numerical method that can deal with non-negative data has been proposed by Y. Brenier et J.D. Benamou through a Fluid Mechanics formulation [START_REF] Benamou | A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem[END_REF].

Fluid Mechanics Formulation

Instead of computing directly the transport, it is possible to consider the geodesic path between the two densities according to the Wasserstein metric (the so-called displacement interpolation [START_REF] Mccann | A convexity principle for interacting gases[END_REF]). For the L 2 ground cost, the geodesic path between the measures with densities ρ 0 and ρ 1 can be shown to have density t → ρ(t, x) where the additional dimension t ∈ [0, 1] parameterizes the path interpolating linearly between the identity map Id d and the Optimal Transport map T :

ρ 0 (x) = ρ(t, T t (x)) |det(∇T t (x))| where T t = (1 -t) Id d +tT.
The geodesic can thus be computed by first obtaining the transport and then making the densities evolve. In [START_REF] Benamou | A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem[END_REF], it has been demonstrated that this geodesic solves the following non-convex problem over the densities ρ(t, x) ∈ R + and a velocity field v(t, x) ∈ R d checking the continuity equation.

Theorem 1 (Benamou-Brenier). In case p = 2 the Wasserstein distance between ρ 0 and ρ 1 is such that:

W (ρ 0 , ρ 1 ) = inf ρ,v Ω 1 0 ρ(t, x)|v(t, x)| 2 dxdt (2.6)
the infimum being taken on ρ, v verifying the non-linear constraint

C 0 = {(ρ, v), ∂ t ρ + div x (ρv) = 0, v, n = 0, ρ(0, •) = ρ 0 , ρ(1, •) = ρ 1 }, (2.7)
with homogeneous Neumann conditions on the velocity field v, through n which is the normal of the domain Ω.

While the estimation of the L 2 Wasserstein distance (2.4) is a static problem over Ω, the introduction of the additional dimension t and the velocity field v leads to a dynamic problem involving kinetic energy (2.6) and continuity equation (2.7) on Ω × [0; 1].

Convexification

Following [START_REF] Benamou | A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem[END_REF] and introducing the change of variable (ρ, v) → (ρ, m), where m is the momentum m = ρv, we obtain a convex optimization problem over the couple (ρ, m):

min (ρ,m)∈C J (ρ, m) := 1 0 Ω J(ρ, m)dxdt (2.8)
which deals the positivity constraint [START_REF] Cardaliaguet | Geodesics for a class of distances in the space of probability measures[END_REF]:

J(ρ, m) =    ||m|| 2 2ρ if ρ > 0 0 if (ρ, m) = (0, 0) +∞ otherwise.
(2.9)

The set of constraint C becomes linear and reads:

C(ρ, m) = {(ρ, m), ∂ t ρ + div x m = 0, m, n = 0, ρ(0, •) = ρ 0 , ρ(1, •) = ρ 1 }. (2.10) Denoting as ι C (x) = 0 if x ∈ C, +∞ otherwise, the convex problem introduced in [19] reads: (ρ * , m * ) = argmin ρ,m J (ρ, m) + ι C (ρ, m). (2.11)
Notice that with this formulation, we are not able to compute the transport map and we only estimate the geodesic ρ(x, t). The existence of minimizers of such a problem in the space of measures has been studied in [START_REF] Villani | Optimal Transport: Old and New[END_REF]. The existence of L 2 minimizers for positive data (ρ 0 , ρ 1 ) has been studied in [START_REF] Guittet | On the time-continuous mass transport problem and its approximation by augmented lagrangian techniques[END_REF]. The existence of L 2 minimizers for nonnegative data (ρ 0 , ρ 1 ) has been shown during the PhD of Romain Hug [START_REF] Hug | Mathematical and numerical analysis of the dynamic optimal transport problem[END_REF], in which the uniqueness of the variable ρ has also been established.

As the problem (2.11) is convex and non-smooth, the Alternating Direction Method of Multipliers (ADMM, see (A.9)) algorithm was considered in [START_REF] Benamou | A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem[END_REF] to compute a global minimizer. Nevertheless, other proximal method that are introduced in Appendix A can be used to compute a global minimizer.

Discretization and Optimization

In order to solve problem (2.11) over the sptaial domain of the image denoted as Ω, the dimension has to be increased with T time steps of artificial time t. The dimension of the variable ρ to estimate is therefore of size |Ω|T .

In [START_REF] Papadakis | Optimal transport with proximal splitting[END_REF], it has been proposed to rely on a staggered grid for the discretization of the continuity equation (2.10), similarly to the discretization of PDE's in incompressible fluid dynamics (see for instance [START_REF] Harlow | Numerical calculation of time-dependent viscous incompressible flow of fluid with free surface[END_REF]). An interpolation operator I has then been considered for computing the cost function (2.9) on the regular grid. The use of a staggered grid is very natural in the context of the discretization of a divergence operator associated to a vector field on R d . The basic idea is to perform an accurate evaluation of every partial derivative at prescribed nodes of a cartesian grid using standard centered finite differences. Hence, the discrete problem (2.11) can be reformulated as:

(ρ * , m * ) = argmin (ρ,m)∈Es J (I(ρ, m)) + ι C (ρ, m), (2.12)
The problem (2.12) consists in minimizing the sum of two non-smooth convex functionals, the first one being composed with a linear operator I. We can thus make use of proximal splitting methods such as Douglas Rachford (A.6) or Primal Dual (A.15) algorithms to solve this problem. As illustrated in the Figure 2.2, faster convergence was observed considering these algorithms instead of the ADMM one. To apply such algorithms that are described in the Appendix, the fundamental point is to be able to compute the proximal operators (defined in relation A.2) of functions J and ι C .

Since the functional J is separable, the proximal operator of J can be computed independently as the proximal operator of J for each point (ρ, m)(t, x) ∈ R × R d with (t, x) ∈ [0; 1] × Ω. As shown in [START_REF] Papadakis | Optimal transport with proximal splitting[END_REF], it simply requires to solve a third order polynomial in each grid point.

Finally, the proximal mapping of ι C is Proj C the orthogonal projector on the convex set C. It requires solving a Poisson equation on the centered grid with prescribed boundary conditions. It can be achieved with Fast Fourier Transform in O(N P log(N P )) operations where N and P are number of spatial and temporal points, see [START_REF] Strang | Linear Algebra and Its Applications[END_REF].

Using this method, one can estimate the geodesic path between densities that vanishes on the spatial domain, as illustrated in Figure 2.1.

t = 0 t = 1/4 t = 1/2 t = 3/4 t = 1

Limitations and Motivations

A first limitation of this approach concerns the computational time. As exhibited in Figure 2.2, as the problem is not strictly convex, the convergence of the cost function J is fast but the convergence of the iterates (ρ k , m k ) is much slower.

(a) Numerical instabilities may also happen when input data are not smooth enough, as illustrated by Figure 2.3. This is nevertheless not a issue for the targeted application where data are images of the ocean representing continuous fluids.

J (ρ k , m k ) (b) ||ρ k -ρ * || 2
ρ 0 ρ(1/6) ρ(1/3) ρ(1/2) ρ(2/3) ρ(5/6) ρ 1 Figure 2.
3 -Image interpolation between Gaspard Monge (ρ 0 ) and Leonid Kantorovitch (ρ 1 ). Numerical instabilities appear (mainly on the boundaries of the head) since the data are not smooth enough.

Next, when solving the Optimal Transport problem between two densities, we know that the mass transfer will follow straight lines. The structures contained in the original data can thus disappear along the optimal path. Such behavior is illustrated in Figure 2.4 that presents the optimal path between images containing two Gaussians. One main objective for image interpolation is to incorporate some physical priors in order to preserve the structures contained in the data along the optimal path. This will be of main interest for geoscience imaging applications. Ocean and atmosphere, observed with satellites, are indeed driven by complex physical laws. The temporal interpolation of such satellite images is an important problem in this community and it should correspond to the underlying dynamics. As interpolation with Optimal Transport involves constant mass transport along straight lines, it is not a satisfactory solution.

ρ(0) = ρ 0 ρ(1/4) ρ(1/2) ρ(3/4) ρ(1) = ρ 1
In order to be able to process oceanographic images that contains "obstacles" such that continent or islands, the topology of the considered spatial model should be carefully taken into account. Some generalizations are proposed in the next chapter to tackle such issues.

Chapter 3

Generalizations for Image Interpolation

In this Chapter, we consider generalized distances based on convex ground costs within the dynamical Optimal Transport formulation. In the aim of interpolating satellite images of the ocean in the presence of coasts, the computation of Optimal Transport on Riemannian manifolds modeling the spatial domain is tackled. Physical regularizations of the transport map are also considered. The numerical scheme proposed in the Chapter 2 can be used with only slight modifications with respect to the L 2 -Wasserstein case.

Generalized Cost Functions

The formulation of the geodesic computation as a convex optimization problem presented in Chapter 2 enables the definition of various metrics obtained by changing the objective function. Transport with congested dynamics [START_REF] Buttazzo | An optimization problem for mass transportation with congested dynamics[END_REF] or unbalanced transport between densities that have different mass [START_REF] Lombardi | Eulerian models and algorithms for unbalanced optimal transport[END_REF][START_REF] Chizat | An Interpolating Distance between Optimal Transport and Fischer-Rao[END_REF] are namely possible. One can also interpolate between different distances. An interpolation between the L 2 -Wasserstein and L 2 distances is proposed in [START_REF] Benamou | Wasserstein optimal mapping between prescribed densities functions[END_REF]. Lastly, an interpolation between L 2 -Wasserstein and H -1 distances is described in [START_REF] Dolbeault | A new class of transport distances between measures[END_REF]. This extension relies in a crucial manner on the convexity of the extended objective function, which enables a theoretical analysis to characterize minimizing geodesics [START_REF] Cardaliaguet | Geodesics for a class of distances in the space of probability measures[END_REF].

Optimal Transport on Riemannian manifolds Many properties of the L 2 -Wasserstein distance can be extended to the setting where the ground cost is the square of the geodesic distance on a Riemannian manifold. This includes in particular the existence and uniqueness of the transport map [START_REF] Mccann | Polar factorization of maps on Riemannian manifolds[END_REF]. Displacement interpolation for transport on manifolds has the same variational characterization as the one introduced in [START_REF] Benamou | Wasserstein optimal mapping between prescribed densities functions[END_REF] for Euclidean transport, see [START_REF] Villani | Optimal Transport: Old and New[END_REF] for a detailed review of Optimal Transport on manifolds.

Displacement interpolation between two measures, each one composed of a single Dirac, amounts to computing a single geodesic curve on the manifold. Discretization and numerical solutions to this problem are numerous. A popular method is the Fast Marching algorithm introduced jointly by [START_REF] Sethian | A fast marching level set method for monotonically advancing fronts[END_REF][START_REF] Tsitsiklis | Efficient algorithms for globally optimal trajectories[END_REF] for isotropic Riemannian metrics (i.e. when the metric at each point is a scalar multiple of the identity) discretized on a rectangular grid. This algorithm has been extended to compute geodesics on 2-D triangular meshes with only acute angles [START_REF] Kimmel | Computing geodesic paths on manifolds[END_REF]. More general discretizations require the use of slower iterative schemes, see for instance [START_REF] Bornemann | Finite-element discretization of static hamiltonjacobi equations based on a local variational principle[END_REF].

Interpolation between pairs of measures in Riemannian manifolds generalizes to barycenters of a family of measures, see [START_REF] Kim | Multi-marginal optimal transport on Riemannian manifolds[END_REF]. The numerical computation of Optimal Transport on manifolds has been less studied. For weighted sums of Diracs, displacement interpolation is achieved by solving a linear program to compute the coupling between the Diracs and then advancing the Diracs with the corresponding weights and constant velocity along the geodesics.

Hence, we have proposed in [START_REF] Papadakis | Optimal transport with proximal splitting[END_REF] to extend the method of [START_REF] Benamou | Wasserstein optimal mapping between prescribed densities functions[END_REF] to solve for the displacement interpolation on a Riemannian manifold. Following [START_REF] Dolbeault | A new class of transport distances between measures[END_REF][START_REF] Cardaliaguet | Geodesics for a class of distances in the space of probability measures[END_REF], we use a generalized cost function that allows one to compute geodesics that interpolate between the L 2 -Wasserstein and the H -1 geodesics . To introduce further flexibility, we have also considered a spatially varying tensor [START_REF] Hug | Multiphysics optimal transportation and image interpolation[END_REF], which corresponds to approximating a transportation problem on a Riemannian manifold.

To that end, the convex cost functional (2.9) is generalized, for β ∈ [0; 1] and A(t, x) a symmetric positive definite tensor, as

J A β (ρ, m) =    Am 2 2ρ β if f > 0, 0 if (m, ρ) = (0, 0), +∞ otherwise. (3.1)
and we now consider the problem:

(ρ * , m * ) = argmin ρ,m 1 0 Ω J A β (ρ, m) dx d t + ι C (ρ, m). (3.2)
The matrix A(t, x) is of size d × d and represents the anisotropic penalization of the displacement energy. The isotropic case of constant weights A = Id d is studied in [START_REF] Dolbeault | A new class of transport distances between measures[END_REF][START_REF] Cardaliaguet | Geodesics for a class of distances in the space of probability measures[END_REF].

The case β = 1 corresponds to the Wasserstein L 2 distance. In a continuous (not discretized) domain, the value of the problem (3.2) for β = 0 is equal to the H -1 Sobolev norm over densities ρ 0 -ρ 1 H -1 , as detailed in [START_REF] Dolbeault | A new class of transport distances between measures[END_REF]. In this case, the induced distance is an Hilbertian norm, and the corresponding geodesic is a linear interpolation of the input measures. Thus, for measures having densities, one obtains ρ t = (1 -t)ρ 0 + tρ 1 .

Interpolation between L 2 -Wasserstein and H -1

With this previous formulation (3.1), it is possible to interpolate between L 2 -Wasserstein and H -1 geodesics. This is illustrated in With this β parameter, smoother interpolation are then obtained. As illustrated in Figure 3.2, it is of interest for interpolating between non-smooth images and avoid numerical instabilities.

β = 0 β = 1/4 β = 1/2 β = 3/4 β = 1 t = 0 t = 1/8 t = 1/4 t = 3/8 t = 1/2 t = 5/8 t = 3/4 t = 7/8 t = 1
β = 1 β = 3/4 ρ 0 ρ(1/6) ρ(1/3) ρ(1/2) ρ(2/3) ρ(5/6) ρ 1

Modelling of obstacle

When β = 1 and the matrices A(t, x) = Id d w(x) are diagonal and constant in time, the solution of (2.12) discretizes the displacement interpolation between the densities (ρ 0 , ρ 1 ) for a ground cost being the squared geodesic distance on a Riemannian manifold.

We exemplify this setting by considering Optimal Transport with obstacles, which corresponds to choosing weights w that are infinity on the obstacle

O ⊂ R × R d , i.e. w(t, x) = 1 + ι O (t, x) ∈ {1, +∞}.
Note that with such definition, the obstacles can be dynamic, i.e. the weight w does not need to be constant in time. Figure 3.3 shows a first example where O is a 2-D (d = 2) static labyrinth map (the walls of the labyrinth being the obstacles are displayed in black). We use a 50 × 50 × 100 discretization grid of the space-time domain [0, 1] 3 and the input measures (ρ 0 , ρ 1 ) are Gaussians with standard deviations equal to 0.04.

For Gaussians with such a small variance, this example shows that the displacement interpolation is located closely to the geodesic path between the centers of the gaussians. Isotropic Riemannian metrics corresponds to matrices A(t, x) = w(x) Id d proportional to the identity at each point, but this extends to arbitrary Riemannian metrics A(t, x) = A(x). The existence of a solution of the Monge-Kantorovitch problem on Riemannian manifolds has been shown in [START_REF] Mccann | Polar factorization of maps on Riemannian manifolds[END_REF]. The equivalence with a fluid mechanism formulation has been studied in [START_REF] Hug | Multiphysics optimal transportation and image interpolation[END_REF]. In the case of a time dependent tensor A(t, x), even if the numerical algorithms are still providing a seemingly sound solution (see Fig- Figure 3.4 shows a more complicated setting that includes a labyrinth with moving walls: a green wall appears at time t = 1/4 and a red one disappears at time 1/2. The difference with respect to the previous example is the fact that w is now time dependent. This simple modification has a strong impact on the displacement interpolation. Indeed, the speed of propagation of the mean of the density is not constant anymore since the density measure is confined in a small area surrounded by the walls for t ∈ [1/4, 1/2].

t = 0 t = 1/9 t = 2/9 t = 1/3 t = 4/9 t = 5/9 t = 2/3 t = 7/9 t = 8/9 t = 1

Anisotropic transport

The introduction of the anisotropy modeling through non diagonal matrices A is now illustrated. With adequate matrices A, it is now possible to define a polarization of the space for modeling anisotropic domain priors. This problem is illustrated in Figure 3.5, which presents the optimal path between a horizontal and a vertical line on a square image. For these experiments, we defined different complex anisotropies, that are represented with the arrows of the first image of each line. If prior knowledge is available on the data, one can build a specific Riemannian manifold and the anisotropic model will therefore be able to simulate rigid or divergence free transports. This modelling is nevertheless not sufficient for more general purposes and we now propose to directly include physical prior into the model and not into the domain.

t = 0 t = 1/9 t = 2/9 t = 1/3 t = 4/9 t = 5/9 t = 2/3 t = 7/9 t = 8/9 t = 1
ρ 0 ρ * (1/8) ρ * (1/4) ρ * (3/8) ρ * (1/2) ρ * (5/8) ρ * (3/4) ρ * (7/8) ρ 1

Optimal Transport with physical priors

With the previous formulation, it is now possible to deal with islands for processing ocean images. Our underlying objective is to incorporate physic priors into the Optimal Transport model. More precisely, for general images, we would like to consider incompressible or rigid transports. Such transports can indeed prevent the object contained within the data from splitting along the computed paths, as for classical transport of Figure 2.4. The physical constraints can be characterized using the velocity, that is not a variable of the problem (2.11) anymore and should be reintroduced. We thus come back to the classic Optimal Transport framework with β = 1 and A = Id d .

Non-convex coupling

In order to have a coupling between the variables (ρ, m) and v, a natural idea is to consider the term ι D (ρ, m, v) = 0, with the set D = {m = ρv}. As D is not convex and ι D non-smooth, we rather consider a differentiable coupling1 :

K(ρ, m, v) = 1 2 Ω 1 0 ||m -ρv|| 2 dx d t.
In [START_REF] Hug | Multiphysics optimal transportation and image interpolation[END_REF], specific kind of transports were promoted regarding velocity priors that depend on the targeted application, as can be done with physical regularization of Optical Flow for fluid image sequences [START_REF] Corpetti | Dense estimation of fluid flows[END_REF][START_REF] Yuan | Discrete orthogonal decomposition and variational fluid flow estimation[END_REF]. A new functional term R(v) is thus introduced. This term can be defined as:

• Divergence-free:

R(v) = ι C v (v), C v = {div x v(t, •) = 0, v(t, •), n = 0, ∀ t ∈ [0; 1]}, • Incompressible penalization: R(v) = 1 0 || div x (v(t, •))|| 2 L 2 (Ω) d t, • Rigid penalization: R(v) = 1 0 ||(∇ x v(t, •) + (∇ x v(t, •)) T )/2|| 2 L 2 (Ω) d t, • Translation penalization: R(v) = 1 0 ||(∇ x v(t, •)|| 2 L 2 (Ω) d t.
By translation penalization, we mean a penalization of the deviation of a velocity field from translations. These terms are convex and the three last ones are differentiable. The existence of minimizers for the translation case has been shown in [START_REF] Hug | Multiphysics optimal transportation and image interpolation[END_REF]. The generalized Optimal Transport model we are interested in now reads:

min ρ,m,v F (ρ, m, v) := J (I(ρ, m)) + ι C (ρ, m) + λK(I(ρ, m), I(v)) + αR(v) (3.3)
where λ, α ≥ 0 respectively weight the coupling between variables and the velocity regularization term. The problem (3.3) is not convex in (ρ, m, v) but it is separately convex in (ρ, m) and in v. As the coupling is differentiable and the non-smooth terms are separable, following [START_REF] Tseng | Convergence of a block coordinate descent method for nondifferentiable minimization[END_REF], we can perform a block coordinate descent with Algorithm (A. [START_REF] Benamou | A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem[END_REF]) and minimize alternatively each convex problem in (ρ, m) and in v to obtain a critical point of the joint problem (3.3). Each subproblem can be solved with Algorithms presented in Appendix A.

Synthetic tests

First of all, we compare in Figure 3.6 the results obtained with the incompressible, translation and rigid penalizations on the two Gaussians example. Contrary to standard transport shown in Figure 2.4, such physical priors prevent the mass from splitting along the computed path. The Gaussians are deformed with the divergence-free prior, but it can be seen that the length of the level lines of the densities are preserved along the path.

It is also important to underline that both translation and rigid penalizations keep the exact shapes of the two Gaussians along time. One can see in Figure 3.7 representing both computed paths, that the rigid penalization really performs a rotation and not a translation, so that the optimal path is no more composed of straight lines. In the example of Figure 3.8 that presents a rotating bar, the rigid penalization (last line) recovers a quasi-rotation, which better preserves the prior physics with respect to pure Optimal Transport (first line). As expected, it can also be observed in Figure 3.9 that the length of the level lines of the estimated density are preserved with the incompressible and rigid penalization approaches. We refer the reader to [START_REF] Bouharguane | Multiphysics optimal transportation and image analysis[END_REF][START_REF] Hug | Multiphysics optimal transportation and image interpolation[END_REF]

Image interpolation in oceanography

In order to study the state of the oceans, snapshot images are produced by operational numerical codes such as the Ocean Circulation Model NEMO (http://www.nemoocean.eu/) . The synthesis process is nevertheless time-consuming and the oceanographers would like to create a few images and then realize a temporal interpolation between these images. This would be of particular interest to visualize and diffuse movies of dynamic structures of the ocean. The main issue comes from the coast that appears in a lot of interesting places, as illustrated in Figure 3.10. It makes useless classical image registration techniques such as optical flow or diffeomorphism estimations that can not deal with such complex domain.

With our formulation, the image domain can be represented as a Riemannian manifold, by taking A(t, x) = w(x) Id 2 . The variable w(x) describes the manifold. It is set to 1 in the ocean and to a huge value in the land in order to restrict the transport into the ocean.

The successive optimal paths between 10 pairs of Sea Surface Height images of size 843 × 516 produced by the NEMO model have been computed with a temporal discretization of 9 steps. The brightest colors correspond to the highest sea height. The image sequence illustrates the Agulhas Current and the creation of vortexes in Cap Point. As shown in Figure 3.10, that presents the computed path between two consecutive images, by adding the proposed divergence-free penalization, we can induce some rotational prior within the transport estimation and recover the creation of vortexes.

Conclusion of Part I

In this first part, we have shown how proximal splitting schemes offer an elegant and unifying framework to describe computational methods to solve the dynamical Optimal Transport with an Eulerian discretization. It allowed us to extend the original method of [START_REF] Benamou | A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem[END_REF] in several directions, most notably the use of staggered grid discretization and the introduction of generalized, spatially variant, cost functions.

We have studied generalized Optimal Transport models which attach a multiphysics model to the images to be registrated. This is of particular interest for image interpolation purposes, where the results obtained using a simple minimization of a kinetic energy under some constraints do not preserve image characteristics along the optimal interpolation path, which is not physical. These generalizations are not limited by the expressions we consider here and others physical terms can also be considered taking into account more complex physics.

Promising results have been obtained on high-resolution oceanographic images. Future works will therefore be dedicated to the modeling of more accurate physical priors.

From the optimization point of view, we also focused on the minimization of nonconvex and non-smooth functionals. Contrary to the block coordinate descent method we considered in (3.3), defining an algorithm without inner loops (on (ρ, m) and v) would be of main interest to obtain a faster computation of minima of our non-convex problems. Recent advances of proximal splitting methods have been made for related problems [START_REF] Attouch | Convergence of descent methods for semialgebraic and tame problems: proximal algorithms, forward-backward splitting, and regularized gauss-seidel methods[END_REF][START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF], but they are still limited to simple functionals and would still require inner loops to deal with our generalized Optimal Transport formulation.

Finally, for the anisotropic Optimal Transport, we plan to study the existence of minimizers when the anisotropy varies with respect to time. This is a challenging problem as, in this case, there are no correspondences with any static Monge problem.

Notice that this method has found interest in the Computer Graphics community. Our original Matlab code has indeed been optimally reimplemented in C++ [START_REF] Bonneel | Faster and specialized c++ implementation of "optimal transport with proximal splitting[END_REF], where 56× speed up has been observed. It is therefore possible to estimate a good geodesic path discretized with 256 timesteps between 256 × 256 grayscale images in 3 minutes. An extension of the proposed model to the interpolation of color images have also been proposed in [START_REF] Fitschen | Dynamic optimal transport with mixed boundary condition for color image processing[END_REF].

Part II

Discrete Optimal Transport

Chapter 4

Relaxation and Regularization of Optimal Transport

An easiest way to discretize and compute numerically Optimal Transport is to consider finite sums of weighted Diracs. In this specific case, the Optimal Transport is a multivalued map between the Diracs locations. The dimension of the transport map then scales with the product of the dimensions of the discretized densities. Despite being numerically intensive for finely discretized distributions, this discrete transport framework has found many applications. It first includes image retrieval with the well-kwown Earth Mover's Distance [START_REF] Rubner | A metric for distributions with applications to image databases[END_REF], and also color transfer between images [START_REF] Pitié | The linear monge-kantorovitch linear colour mapping for example-based colour transfer[END_REF], shape retrieval [START_REF] Rabin | Geodesic shape retrieval via optimal mass transport[END_REF] or surface reconstruction [START_REF] De Goes | An optimal transport approach to robust reconstruction and simplification of 2D shapes[END_REF] and interpolation for computer graphics [START_REF] Bonneel | Displacement interpolation using Lagrangian mass transport[END_REF].

As already seen in the previous part, the Optimal Transport map between complicated densities is usually irregular, which is a real drawback for interpolating between densities or for transfer purposes. In this chapter, we first recall the discrete setting of Optimal Transport and next present some regularization approaches that will be used for image processing purposes in Chapters 5 and 6.

Discrete Formulation

Monge's original formulation of the Optimal Transport problem aims at minimizing the cost for transporting a distribution µ onto another distribution ν using a map T min

T X c(x, T (x))dµ(x), where T #µ = ν. (4.1) Here, µ, ν are measures in R d , T : R d → R d is a µ-measurable function, c : R d ×R d → R + is a µ ⊗ ν Y -measurable function,
and # is the push forward operator. We now consider the discrete formulation of the optimal mass transportation problem between a pair of histograms µ and ν. In the most general setting, histograms can be viewed as weighted sum of dirac distributions

µ = N X i=1 µ i δ X i and ν = N Y j=1 ν j δ Y j ,
where δ X is the Dirac measure at location X ∈ R d and

X = {X i ∈ R n } N X i=1 and Y = {Y j ∈ R n } N Y
j=1 are point-clouds in a feature space R d (for instance position, color, gabor filter coefficients, etc). These points are often referred to as "bins" when considering a regular grid. Without loss of generality, we assume from now that the two histograms have positive weights µ i , ν j > 0, the same number of bins (N X = N Y = N ) and the same total mass i µ i = j ν j .

Monge-Kantorovitch Formulation

Given a fixed cost matrix C = (C i,j ) i,j measuring a distance between locations X i and Y j , we can consider the Kantorovitch formulation of the Optimal Transport problem, which corresponds to estimating the Optimal Transport plan P * between histograms. The Monge-Kantorovitch distance MK reads [START_REF] Kantorovitch | On the translocation of masses[END_REF] MK(µ, ν) = P , C = min

P ∈S(µ,ν) P, C = N i=1 N j=1 P i,j C i,j , (4.2)
where the set of admissible transport matrices S(µ, ν) is the set of measures defined on the product space of considered histograms which marginals are given by µ and ν:

S(µ, ν) := {P ∈ R N ×N + , P 1 N = µ and P T 1 N = ν}. (4.3)
As illustrated in Figure 4.1, the value P i,j of matrix P then corresponds to the proportion of the mass µ i of cluster X i that is transferred to the cluster Y j . If C i,j = X i -Y j p , the value of the optimization problem (4.2) is called the L p -Wasserstein distance (to the power p) and is denoted W p (µ, ν) p = MK(µ, ν). It can be shown that W p defines a distance on the set of distributions that have moments of order p. Problem (4.2) can be reformulated as a linear program that scales with the product of the dimensions of the consider histogram. The computation of the transport cost is therefore limited in practice to histograms discretized with a low number of bins.

Linear programming methods (interior-point method [START_REF] Nesterov | Interior Point Polynomial Methods in Convex Programming: Theory and Algorithms[END_REF], simplex algorithm [START_REF] Dantzig | Linear Programming and Extensions[END_REF]) or first-order methods (Proximal Algorithms presented in Appendix A, Frank-Wolfe algorithm, as proposed in [START_REF] Zaslavskiy | A path following algorithm for the graph matching problem[END_REF]) can be applied to estimate the transport matrix P . Faster schemes exist for specific cost functions C X,Y , such as for instance convex cost of the distance on the line (where it boils down to sorting the positions as d = 1) and the circle [START_REF] Delon | Fast transport optimization for Monge costs on the circle[END_REF], concave costs on the line [START_REF] Delon | Local matching indicators for transport problems with concave costs[END_REF], the 1 distance [START_REF] Ling | An efficient earth mover's distance algorithm for robust histogram comparison[END_REF].

In case d > 1, various approximations of the transportation distance have been proposed using Kantorovitch-Rubinstein discrepency [START_REF] Shirdhonkar | Approximate earth mover's distance in linear time[END_REF] or iterative explicit 1D computations [START_REF] Bonneel | Sliced and radon wasserstein barycenters of measures[END_REF]. The computation can also be accelerated with multi-scale clustering [START_REF] Mérigot | A multiscale approach to optimal transport[END_REF][START_REF] Schmitzer | A hierarchical approach to optimal transport[END_REF][START_REF] Oberman | An efficient linear programming method for Optimal Transportation[END_REF], using the structure of the cost function [START_REF] Schmitzer | A sparse algorithm for dense optimal transport[END_REF] or considering Voronoi diagrams [START_REF] Aurenhammer | Minkowski-type theorems and least-squares partitioning[END_REF][START_REF] Levy | A numerical algorithm for $l_2$ semi-discrete optimal transport in 3d[END_REF].

Remark 1 (Transport between clouds of points). In the case where the measures are discrete, have the same number of points, and all points have the same mass, µ and ν are cloud of points and the transport T between them is a one-to-one assignment, which corresponds to optimizing over the set of permutation matrix P in the Kantorovitch formulation (4.2). Standard combinatorial methods such as the Hungarian Algorithm [START_REF] Kuhn | The Hungarian method of solving the assignment problem[END_REF] or the Auction Algorithm [START_REF] Bertsekas | The auction algorithm: A distributed relaxation method for the assignment problem[END_REF] can thus be used to compute the transport matrix. Notice that the set of permutation matrices is not convex, so one can rely on its convex hull, the set of bi-stochastic matrices. It is shown that the relaxation over this set is tight, meaning that there exists a solution of (4.2) which is a binary matrix, hence being also a solution of the original non-convex problem over permutation matrices, see [START_REF] Villani | Topics in Optimal Transportation[END_REF].

Dual formulation

An alternative for solving the Optimal Transport problem (4.2) consists in looking at its dual formulation:

MK(µ, ν) = max u∈R N , v∈R N s.t. u i +v j ≤C i,j , ∀ i,j µ T u + ν T v, (4.4)
which comes from the Legendre-Frenchel transform of the primal problem:

MK * (u, v) = ι u i +v j ≤C i,j (u, v) (4.5)
where ι K stands for the characteristic function of the convex set K. In this case the number of unknowns boils down to 2N but the number of linear constraints is N 2 . Hence, solving the dual problem is as hard as solving the primal one. This dual formulation presents some advantages for particular image processing applications. It can for instance be considered to convexify vecto-valued labelling problems [START_REF] Strekalovskiy | Convex relaxation of vectorial problems with coupled regularization[END_REF].

Relaxed and Regularized Transport

For applications such as color transfer, the regularization of the transport map is a necessary ingredient to provide visually plausible images. To that end, a simple idea consists in adding a regularization penalty to the Optimal Transport energy (4.2). However, it leads to difficult non-convex problems, that have not yet been solved in a satisfying manner either theoretically or numerically.

Graph regularization and matching Based on concepts developed in manifold learning [START_REF] Tenenbaum | A global geometric framework for nonlinear dimensionality reduction[END_REF], regularizations built on top of a graph structure connecting neighboring points in the input density are classically used in imaging applications [START_REF] Elmoataz | Nonlocal discrete regularization on weighted graphs: A framework for image and manifold processing[END_REF]. With graphs, one can design regularizations that are adapted to the geometry of the input density, that often has a manifold-like structure.

This idea of graph-based regularization of Optimal Transport can be interpreted as a soft version of the graph matching problem, which is at the heart of many computer vision tasks, see [START_REF] Belongie | Shape matching and object recognition using shape contexts[END_REF][START_REF] Zheng | Robust point matching for nonrigid shapes by preserving local neighborhood structures[END_REF]. Graph matching is a quadratic assignment problem, known to be NP-hard. Several convex approximations have been proposed, including for instance linear programming [START_REF] Almohamad | A linear programming approach for the weighted graph matching problem[END_REF] and SDP programming [START_REF] Schellewald | Evaluation of a convex relaxation to a quadratic assignment matching approach for relational object views[END_REF].

Transport relaxation In the recent work of [START_REF] Louet | A sharp inequality for transport maps in via approximation[END_REF], it has been shown that in 1-D, no regularization is possible if one maintains a one-to-one assignment between the two densities. This is a first motivation for introducing a relaxed transport which is not a bijection between the densities. Another (more practical) motivation is that relaxation is crucial to solve imaging problems such as color transfer. Indeed, the color distributions of natural images are multi-modals. An ideal color transfer should match the modes together. This cannot be achieved by classical Optimal Transport because these modes often do not have the same mass.

A typical example is for two images with strong foreground and background dominant colors (thus having bi-modal densities) but where the proportion of pixels in foreground and background are not the same. Such simple examples cannot be handled properly with Optimal Transport. Monitoring the variation of mass between the matched densities requires an appropriate relaxation of the mass conservation constraint. Mass conservation relaxation is related to the relaxation of the bijectivity constraint in graph matching, for which a convex formulation is proposed in [START_REF] Zaslavskiy | Many-to-many graph matching: a continuous relaxation approach[END_REF].

In the previous section, we introduced the Monge-Kantorovitch formulation for the computation of the Optimal Transport between two distributions as the minimization of the energy (4.2). In this section, we modify this energy in order to obtain a regular Optimal Transport mapping. We will also consider the recent entropic regularization of OT proposed in [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF] that presents different and interesting properties, namely in terms of computational cost.

Relaxed Transport

In many applications in imaging, strict mass conservation should be avoided. As a consequence, it is not desirable to impose a one-to-one mapping between the points in µ and ν. The relaxation proposed in [START_REF] Ferradans | Regularized discrete optimal transport[END_REF] allows each point of µ to be transported to multiple points of ν and vice versa. For this purpose, the Optimal Transport problem (4.2) is modified as: min

P ∈Sκ(µ,ν) C X,Y , P (4.6)
where 4 are the parameters of the method. They describe the minimum (k X , k Y ) and maximum (K X , K Y ) proportion of each bin that can be transported. To impose the total amount of mass M transported between the densities, it is necessary to further set the constraint 1 T P 1 = M , where M > 0 is a parameter, basically taken as M = 1 for transporting all the mass from µ to ν. To ensure that S κ is non empty, the condition max(k X , k Y ) ≤ M ≤ min(K X , K Y ) has to be fulfilled. Notice that problem (4.6) can still be solved using standard linear programming algorithms.

S κ (µ, ν) = P ∈ [0, 1] N ×N \ k X µ ≤ P 1 ≤ K X µ, k Y ν ≤ P T 1 ≤ K Y ν, 1 T P 1 = M . (4.7) In this new problem, κ = (k X , K X , k Y , K Y ) ∈ (R + )

Case of cloud of points

In the special case where µ and ν are cloud of points, one can consider µ = ν = 1.

The following proposition shows that an optimal P is binary when the parameters κ are integers. Such a binary P can be interpreted as a set of pairwise assignments between the points in X and Y . Note that this is not true in general when the parameters κ are not integers. 5 , there exists a solution P of (4.6) which is binary, i.e. P ∈ {0, 1} N ×N . This property is due to the fact that binary matrices P ∈ S κ (µ, ν) are totally unimodular and the polytope S κ (µ, ν) has integer vertices for integer values κ. Since there is always a solution of the linear program (4.6) which is a vertex of S κ (µ, ν), it has coefficients in {0, 1} [START_REF] Schrijver | Theory of linear and integer programming[END_REF].

Proposition 1. For (k X , K X , k Y , K Y , M ) ∈ (N * )
It is illustrated in Figure 4.2. Given a set of points X (in blue) and Y (in red), the optimal P is computed solving (4.6) for different values of κ. For each values of κ, a line between X i and Y j is drawn if the value of the associated optimal P i,j > 0.1, solid if P i,j = 1, and dashed otherwise. For non integer values of K X , K Y , the mappings P i,j are in [0, 1] while for integer values, P i,j ∈ {0, 1}. As the values of K X , K Y are increased (Figure 4.2, right), the points in X tend to be mapped to the closer points in Y .

Mean Optimal Transport Map

As stated before, when general distributions are concerned, optimal matrices P minimizing (4.6) are in general non binary and their non zero entries do not define one-to-one maps between the points of X and Y . As illustrated in Figure 4.3, it is however possible to define a map T from X to Y by mapping each point X i to a weighted barycenter of its neighbors in Y as defined by P . This corresponds to defining mean that the sum of links connecting each red point to blue points equals 1, while the sum of links connecting each blue point to red points is in [0; 2]. In this case, as all the parameters are integers, assignments are obtained (solid lines) and each red point is assigned to exactly one blue point, whereas each blue point can be connected to 0, 1 or 2 red points.

T (X i ) = N j=1 P i,j Y j N j=1 P i,j
which in vectorial form can be expressed as T (X i ) = ((diag(P 1)) -1 P Y ) i , where the operator diag(v) creates a diagonal matrix in R N ×N with the vector v ∈ R N on the diagonal. To insure that the map is well defined, we impose that k X > 0. Note that it is possible to define a map from Y to X by replacing P by P T in the previous formula and exchanging the roles of X and Y .

Discrete Regularized Transport

So far, a transport problem where the mass conservation constraint is relaxed has been presented. The second step is to define its regularization. A classic way of imposing regularity on a mapping V : R d → R d is to measure the amplitude of its derivatives. Two examples for continuous functions are the quadratic Tikhonov regularizations such as the Sobolev semi-norm ∇V 2 , and the anisotropic total variation semi-norm ∇V 1 regularization. However, the differential operator ∇ cannot be applied directly to a point clouds X i due to the lack of neighborhood definition. To extend the definition of the gradient operator, a graph structure on each point cloud has to be designed.

In the current setting, the goal is to regularize the discrete map T introduced in (4.8), which is only defined at the location of the points as

X i → Ṽi = X i -T (X i ) = X i -diag(P 1) -1 (P Y ) i .
To avoid the normalization diag(P 1) (which typically leads to non-convex optimization problems), the regularity can be imposed on the map

X i → V i = diag(P 1)X i -(P Y ) i .
This switch also has the advantage of imposing a stronger regularization in regions with large weights and to further regularize the variations of the weights P 1 ∈ R N . These are actually interesting features to reduce artifacts for imaging applications.

Gradient on Graphs A natural way to define a gradient on a point cloud X is by using the gradient on a weighted graph G X = (X, E X , W X ) where E X ⊂ {1, . . . , N } 2 is the set of edges and W X is the set of weights, W X = (w i,j

) N i,j=1 : {1, . . . , N } 2 → R + , satisfying w i,j = 0 if (i, j) / ∈ E X .
The edges of this graph are defined depending on the application. A typical example is the n-nearest neighbor graph, where every vertex X i is connected to X j if X j is one of the n-closest points to X i in X, creating the edge (i, j) ∈ E X , with a weight w i,j . Because the edges are directed, the adjacency matrix is not symmetric.

The gradient operator on G X is defined as G X : R N ×d → R P ×d , where P = E X is the number of edges and where, for each

V = (V i ) N i=1 ∈ R d , G X V = (w i,j (V i -V j )) (i,j)∈E X ∈ R P ×d .
A classic choice for the weights to ensure consistency with the directional derivative is w i,j = X i -X j -1 , see for instance [START_REF] Gilboa | Nonlocal linear image regularization and supervised segmentation[END_REF].

Regularity Term

The regularity of a transport map V ∈ R N ×d is measured according to some norm of G X V , that we choose here for simplicity to be the following

J p,q (G X V ) = (i,j)∈Ex w i,j (V i -V j ) q p ,
where . q is the q norm in R d .

The case (p, q) = (1, 1) is the graph anisotropic total variation, (p, q) = (2, 2) is the graph Sobolev semi-norm, and (p, q) = (1, 2) is the graph isotropic total variation, see for instance [START_REF] Elmoataz | Nonlocal discrete regularization on weighted graphs: A framework for image and manifold processing[END_REF] for applications of these functionals to imaging problems such as image segmentation and regularization. In practice, we observe that the optimal matrix P is somehow sparser with the total variation prior, i.e. the resulting weights are more concentrated on a few entries.

Symmetric Regular Optimal Transport Formulation

Given two histograms µ and ν, a relaxed and regularized Optimal Transport mapping between them can be computed. In order to make the mapping regular with respect to both histograms, the displacement fields to be regularized are: T X,Y (P ) = diag(P 1)X -P Y and T Y,X (P * ) = diag(P 1)Y -P T X.

In order to obtain a partial matching that is regular according to X and Y , two graphs G X and G Y are created (as described Page 47) and the corresponding gradient operators are denoted G X ∈ R n X ×N and G Y ∈ R n Y ×N , where n X and n Y are the number of edges in the respective graphs. The symmetric regularized discrete Optimal Transport energy proposed in [START_REF] Ferradans | Regularized discrete optimal transport[END_REF] is defined as:

min P ∈Sκ P, C X,Y + λ X J p,q (G X T X,Y (P )) + λ Y J p,q (G Y T Y,X (P )), (4.9)
where (λ X , λ Y ) ∈ (R + ) 2 controls the desired amount of regularity. The case κ = (1, 1, 1, 1) and (λ X , λ Y ) = (0, 0) corresponds to the usual Optimal Transport defined in (4.2), and (λ X , λ Y ) = (0, 0) corresponds to the un-regularized formulation (4.6).

By setting the parameters q = p = 1 to obtain the anisotropic Total Variation (TV) norm, problem (4.9) can again be re-written as a linear program by introducing the auxiliary variables U

X ∈ R n X ×d and n Y ∈ R P Y ×d min P,U X ,U Y C X,Y , P + λ X U X , 1 + λ Y U Y , 1 subject to      -U X ≤ G X (P Y -diag(P 1)X) ≤ U X , -U Y ≤ G Y (P T X -diag(P T 1)Y ) ≤ U Y , P ∈ S κ .
(4.10)

In Figure 4.4, the influence of the parameters κ and (λ X , λ Y ), from equation (4.9) are illustrated on a synthetic example. For λ X = λ Y = 0 one obtains the relaxed symmetric Optimal Transport solution, where the transport maps the points in X to the closest point on Y , and vice versa. As we increase the values of λ X and λ Y to 0.001, we can see how the regularization affects the mapping. Let us analyze J p,q (G X T X,Y (P )) = G X diag(P 1)X -G X P Y 2 , for instance. The term G X diag(P 1)X is measuring the regularity of the weights diag(P 1) on X and the consequence is that for λ X = λ Y = 0.001 there are plenty of connections with low weight (there are few solid lines), while for λ X = λ Y = 0 there are several mappings with P i,j = 1 (solid lines). So, the regularization promotes a spreading of the matchings.

The minimum of J p,q (G X T X,Y (P )) is reached when G X diag(P 1)X = G X P Y , that is, when the graph structure of X has the same shape as the graph structure of P Y , which both can be observed in the last column and row. For high values of λ X = λ Y the matchings tend to link the clusters by their shape, that is, the big cluster on X with the big cluster of Y , and similarly for the small clusters (note that the links with higher value are between the small clusters).

We will show in the next chapter how such framework can be efficiently applied to color transfer. In the meanwhile, we finally present the entropic regularization of Optimal Transport that will be mainly used for segmentation purposes. The introduction of regularization (λ X = λ Y = 0.001) spreads the connections among the clusters, while maintaining the cluster-to-cluster matching. For a high value of λ X = λ Y = 10, the regularization tends to match the clusters with similar shape with each other, where the shape is defined by the graph structure. The graphs G X and G Y are represented with the nodes on blue and red respectively, and the edges as solid lines.

λ X = λ Y = 0 λ X = λ Y = 0.001 λ X = λ Y = 10 Graphs

Approximate MK cost using Sinkhorn distances

The entropy smoothing of Optimal Transport proposed in [START_REF] Wilson | The Use of Entropy Maximising Models, in the Theory of Trip Distribution, Mode Split and Route Split[END_REF][START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF] reads

MK λ (µ, ν) = P λ , C = min P ∈S(µ,ν) P, C - 1 λ h(P ) , (4.11)
where the entropy of the matrix P is defined as h(P ) :=i,j P i,j log P i,j . As shown in [START_REF] Rangarajan | A convergence proof for the softassign quadratic assignment algorithm[END_REF], it acts as a barrier function that ensures the positivity of P i,j while making the problem strongly convex. The transport cost MK λ (µ, ν) is called the Sinkhorn distance.

Denoting as K = e -λC , it can be shown [START_REF] Sinkhorn | Diagonal equivalence to matrices with prescribed row and column sums[END_REF][START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF] that the optimal transport map can be written as P λ = diag(u)Kdiag(v) and can be estimated with the following fixed point iteration algorithm on u and v:

u i = µ i j K i,j v j , v j = ν j i K i,j u i (4.12)
One should note that such type of normalization, i.e. x → x x 1 is the projection on the 1 ball (or equivalently probability simplex) using the Kullback-Leibler divergence [START_REF] Benamou | Iterative bregman projections for regularized transportation problems[END_REF]. Indeed, another way to express the regularization term is:

MK λ (µ, ν) = min P ∈S(µ,ν) P, C + 1 λ KL(P 1).
When regular grids X i and Y j are involved and the cost C i,j ∝ X i -Y j , the previous algorithm (4.12) corresponds to convolutions with the kernel K i,j = e -λC i,j . For C i,j = ||X i -Y j || 2 , one recovers convolutions with Gaussian Kernels that can be implemented separately in each dimension thus leading to very fast algorithms [START_REF] Solomon | Convolutional wasserstein distances: Efficient optimal transportation on geometric domains[END_REF].

In Figure 4.5, the matrices P of transport flows computed with the different regularization methods are illustrated. The null values P i,j are displayed in white color. These matrices describe how the mass from the bins of the input histogram (in rows) are sent to the bins of the target histogram (in columns). As a bin µ i corresponds to the line i of the matrix, the variance of the transport map is larger when a lot of columns j are active (i.e. P i,j > 0) and correspond to very different bins Y j . Even if both regularizations have very different origins, we can see on this illustration that the influence over the transport matrix P is quite similar, as the mass of each bin of one histogram is spread to several bins on the other histogram.

The parameter λ of the entropy regularization [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF] directly measures this dispersion. In the asymptotic case λ → 0, the transport matrix is the one with the minimum entropy: P = µν T . In the case λ → ∞, one recovers the standard transport matrix (Figure 4.5 (c)), which is the one of maximum entropy. In practice the convergence of the algorithm (4.12) that estimates the 2N variables is faster for small values of λ. A trade-off between computational cost and quality of the approximate transport has to be chosen.

On the other hand, even if it leads to slower algorithms (since a linear program with N 2 variables has to be solved), the explicit regularization of the transport plan [START_REF] Ferradans | Regularized discrete optimal transport[END_REF] has a physical meaning which presents some advantages as will be seen later.

Dual formulation As shown in [START_REF] Cuturi | Fast computation of wasserstein barycenters[END_REF], the Fenchel transform of (4.11) reads

MK * λ (u, v) = 1 λ i,j e -λ(C i,j -u i -v j )-1 . (4.13)
Note that the indicator function for condition C i,j -αµ i -βν j of the standard transform (4.5) has been smoothed by an exponential function. As noticed in [START_REF] Cuturi | A smoothed dual approach for variational wasserstein problems[END_REF], such function is smooth and its gradient can be computed explicitly, which is a main advantage for incorporating MK terms in more general functionals to be minimized.

More generally, generalized constraints can be naturally integrated with the dual formulation of the entropic regularization. A simple example concerns unnormalized histograms µ and ν. In order to impose them to have a total mass of M , the regularized Optimal Transport function has been extended in [START_REF] Rabin | Convex color image segmentation with optimal transport distances[END_REF] as follows:

MK λ,M (µ, ν) = min P ∈S(µ,ν) P, C - M λ h P M + ι H M (µ, ν), (4.14)
where

H M := {µ, ν ∈ R N + , µ, 1 N = ν, 1 N = M }. (4.15)
The Fenchel transform of the regularized transport cost (4.14) reads:

MK * λ,M (u, v) = M λ log i,j e -λ(C i,j -u i -v j ) , (4.16)
which directly encodes the mass constraint.

In the next chapters, we will see how these different regularizations can be used to deal with image processing problems.

Chapter 5 Application to Color Transfer

Many image processing applications require the modification or the prescription of some characteristics (colors, frequencies or wavelet coefficients) of a given image, while preserving other features. Statistics to be prescribed may come from prior knowledge, or more generally, are learned from an example. In such a case, another image is selected from a database to define a template. Such a framework arises for image enhancement, inpainting, colorization of grayscale or infrared images, tone mapping, color grading or color transfer. In this paper, we will focus on this last application through histogram transfer between images.

Color transfer consists in modifying an image to match the color palette of another one, while preserving its geometry. In the literature, the different interpretations and definitions of color palettes have led to various algorithms. In the following, we only consider unsupervised approaches. Since the seminal work of [START_REF] Reinhard | Color transfer between images[END_REF], methods have been designed to transfer some simple color statistics (i.e. the mean and standard deviation) in any color space [START_REF] Xiao | Color transfer in correlated color space[END_REF]. These works share some common features with texture synthesis by example (see in particular [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF]). More general approaches match the complete histogram of features from two images. When considering grayscale images, the problem is known as 1-D histogram specification. This framework has been extended for color histograms, using for instance 1-D ∆E-color index [START_REF] Morovic | Accurate 3d image colour histogram transformation[END_REF].

Histogram matching via optimal transportation

There exist strong links between histogram specification and the Monge-Kantorovich Optimal Transport problem. Histogram transfer between grayscale images indeed corresponds to the application of the 1-D Optimal Transport plan between the gray levels of the pixels of the images, as noted by [START_REF] Delon | Midway image equalization[END_REF] for the case of histogram equalization. Some approaches to find fast approximate solution of OT for color images were investigated in [START_REF] Pitié | The linear monge-kantorovitch linear colour mapping for example-based colour transfer[END_REF][START_REF] Rabin | Wasserstein barycenter and its application to texture mixing[END_REF].

Spatial information As illustrated in Figure 5.1 (c), the exact transfer of color palette is generally not satisfying for practical applications in image processing [START_REF] Pouli | Progressive color transfer for images of arbitrary dynamic range[END_REF]. The color distributions to be matched may have very different shapes, so that outliers generally appear in the processed image. Moreover, as the process is performed in the color space, it may not transfer coherent colors to neighboring pixels. It results in undesirable artifacts, such as JPEG compression blocks, enhanced noise, saturation or contrast inversion that should be corrected with specific post-processings [START_REF] Pitié | The linear monge-kantorovitch linear colour mapping for example-based colour transfer[END_REF][START_REF] Rabin | Removing artefacts from color and contrast modifications[END_REF][START_REF] Su | Corruptive artifacts suppression for example-based color transfer[END_REF]. As a consequence, various models have been designed to incorporate some regularity priors on the image domain, such as Total Variation [START_REF] Nikolova | Exact histogram specification for digital images using a variational approach[END_REF]. Color transfer may be formalized [START_REF] Papadakis | A variational model for histogram transfer of color images[END_REF] as a variational problem in the image domain, in order to directly incorporate a spatial regularization of colors. While spatial regularization suppresses small artifacts due to exact histogram specification (see Figure 5.1 (d)), it cannot handle strong artifacts due to an irregular OT map [START_REF] Rabin | Wasserstein regularization of imaging problem[END_REF] and deal with potential very different modes (i.e. the dominant colors) between each distribution. For instance, one would like to only transfer a limited quantity of purple color into the flower example of Figure 5.1.

Hence, we investigate in this chapter the use of relaxed and regularized Optimal Transport for proposing models robust to color outliers that directly include spatial information in the image domain.

Convex and regularized Optimal Transport model

This section shows how the relaxed and regularized Optimal Transport formulation introduced in Chapter 4 can be applied to color transfer, and how the regularization and the relaxation improve the results obtained by previous methods.

The color transfer problem consists in modifying a source image I X so that its colors match the colors of a target image I Y . We denote as I X : Ω ⊂ Z 2 → Σ ⊂ R 3 , where Ω is the regular pixel grid and Σ is the quantized 3D RGB color space. The problem is to find a new image I Z whose geometry is as close as possible to the source image I X and whose color distribution is close to the one of the exemplar image I Y . The color histogram of such an image I X can be estimated using the empirical distribution µ X . The goal of color transfer algorithms is to compute a transformation T such that for all pixel x ∈ Ω, I Z (x) = T (I X (x)), where the new empirical distribution µ Z is close (or equal) to µ Y . Problem Formulation We consider the bipartite weighted matching problem between two histograms. Let µ and ν be two discrete and normalized distributions of features X ∈ R N X ×d and Y ∈ R N Y ×d , respectively. Taking N X = N Y = N , their corresponding weights are µ ∈ R N and ν ∈ R N , so that: µ = i µ i δ X i and ν = j ν j δ Y j . We look for a mapping T : R d → R d between X and Y , referred to as transfer map in the following. In order to transport µ towards ν, we define T #µ = i µ i δ T (X i ) . For color transfer, the mapping should have the following properties: (i) T #µ should be close to or match the target distribution ν, (ii) T should be regular in some sense.

Adaptive OT

Relaxing the mass conservation constraint is crucial in order to better match the modes (i.e. the dominant colors) of each distribution. Regularizing the transport is also important to reduce colorization artifacts. Inspired from functional (4.9), the regularized model proposed in [START_REF] Rabin | Adaptive color transfer with relaxed optimal transport[END_REF] aims at tackling the aforementioned issues by relaxing the matching constraint towards the target distribution ν on the matrix P . This is done by introducing the set of acceptable transfer maps P κ (µ, ν):

P κ (µ, ν) = {P ∈ [0, 1] N ×N \ P 1 N = µ, P T 1 N = κ ⊗ ν} where κ ∈ R N , κ ≥ 0, κ ⊗ ν = (κ j ν j ) j and s.t. κ, ν ≥ 1.
Observe here that the second constraint involves a set of local relaxation parameters κ = (κ j ). With respect to the previous admissible set (4.7), the relaxation with these parameters allow the model to increase or decrease the capacity value of the corresponding feature Y j of the target histogram µ Y ; the maximum amount of mass from µ that Y j can receive is now bounded by κ j ν j . If κ = 1 N Y , we get the previous Optimal Transport problem. The last constraint κ, ν ≥ 1 prevents the set of acceptable transport maps P κ (µ, ν) from being empty. With this relaxed parameterization, a model that locally adapt itself to the data can be defined.

Convex formulation When considering color transfer, we would like the transport to be piecewise constant, since we do not want to penalize color shift. To that end, as in (4.8), we consider the regularization of the average transport displacement T (X i ) = (D µ P Y -X) i , where D µ = diag(µ -1 ). As in the previous chapter, this regularization is done on the gradient G X of the average transport T , that is computed on the graph created from the nodes X i . Moreover, in order to control that the relaxed and regularized transported color palette is close to the target, an additional term κ -1 N is added to check that the relaxation is tight, i.e. that the final transported density will be close enough to the proportion of colors in ν.

Introducing the regularization parameter λ ≥ 0 and the fidelity one ρ ≥ 0, the model derived from (4.9) is the following convex optimization problem with linear constraints:

{P , κ } ∈ argmin P ∈Pκ(µ,ν) κ∈R N ,κ≥0, κ, ν ≥1 C XY , P + λ G X (D µ P Y -X)) 1 + ρ κ -1 1 .
(5.1)

This energy can be minimized with a linear solver by introducing auxiliary variables to deal with the L 1 norms, as detailed in (4.10). This formulation generalizes the Optimal Transport model proposed in [START_REF] Ferradans | Regularized discrete optimal transport[END_REF] on cloud of points, to handle distributions. Observe that the transport is here automatically relaxed, since we do not require the user to tune the capacity vector κ.

Spatio-color clustering

To make the optimization problem (5.1) tractable for histograms obtained from large scale images, a joint clustering is performed on the pixel grid and on the color space (d = 5). The images I X and I Y are quantized to generate X and Y with N clusters. We denote X i = (x i , c i ) ∈ R 5 , to specify the spatial component (x i ∈ Ω) and the color component (c i ∈ Σ). These values (x i , c i ) represent the mean spatial and color values of the set of pixels in I X (and equivalently for I Y ) assigned to cluster X i . As illustrated in Figure 5.3 (a-b), this assignment is computed using the fast super-pixels method [START_REF] Achanta | SLIC superpixels compared to state-of-the-art superpixel methods[END_REF]. Note that each cluster X i has an assigned group of pixels {N i ⊂ Ω × Σ} i that can have a different amount |N i | of pixels. We define the cluster weight as

µ i = |N i |/|Ω|.
Thus, we consider from now the measure µ : Graph With this clustering involving spatial information of color location, color consistency in the image domain can be ensured through graph regularization. A graph is created from the nodes X i that are assumed to be linked to P nodes X j . The gradient

X = (x, c) → i µ i δ X i (X) = i µ i δ x i (x)δ c i (c
G X V ∈ R N ×P ×d of a field V = {V l i } ∈ R N ×d on G X is computed at point X i as (G X V ) i = (w ij (V i -V j )) j) ∈ R P ×d ,
where the weight w ij between cluster X i and X j is non null if the nodes are linked in the graph. In this case, the weight relies on the clusters' distances:

w ij ∝ exp(-||X i -Y j ||).
Image Synthesis As we work at a super-pixel scale to speed-up the Optimal Transport computation, the last step is to synthesize a new full resolution image I Z from the source image I X using the new color palette T #µ. Following [START_REF] Tai | Local color transfer via probabilistic segmentation by expectation-maximization[END_REF], a maximum likelihood estimation can be used to incorporate geometrical information from the source image I X into the synthesis process. For each pixel X = (x, I X (x)) of the source image, the idea is to use the likelihood to all clusters {X i } i through the weights {ω i } i in order to compute a linear combination of the estimated transferred colors T (X i ). Hence, we finally get for each pixel x ∈ Ω: (x,

I Z (x)) = Z i (x) = 1 W ( X) i w i (x)T (X i ), with w i (•) = exp -1 2 • -X i 2
, and the normalizing factor W ( X) = i w i ( X).

Color transfer experiments

A first example is presented in

Photo credit: Nicolas Le Dilhuit

A comparison with state-of-the art methods [START_REF] Pitié | The linear monge-kantorovitch linear colour mapping for example-based colour transfer[END_REF], [START_REF] Papadakis | A variational model for histogram transfer of color images[END_REF][START_REF] Ferradans | Regularized discrete optimal transport[END_REF] is given in Figure 5.4. Notice that the simplified optimal transport of colors with Principal Component Analysis of [START_REF] Pitié | The linear monge-kantorovitch linear colour mapping for example-based colour transfer[END_REF] has been improved with a post-processing [START_REF] Rabin | Removing artefacts from color and contrast modifications[END_REF]. Thanks to the super-pixel segmentation, the process only requires a few seconds whereas the methods of the literature last a few minutes. Contrary to [START_REF] Pitié | The linear monge-kantorovitch linear colour mapping for example-based colour transfer[END_REF][START_REF] Papadakis | A variational model for histogram transfer of color images[END_REF], the target color palette is not perfectly matched, since the model adapts itself to the source image data. It leads to a transfer visually more plausible, as illustrated in the flower example of Figure 5.5 where the correct proportion of purple has been selected. With respect to [START_REF] Ferradans | Regularized discrete optimal transport[END_REF], the improvement of the spatial regularization involved by the super-pixel clustering removes artifacts that appear for instance in the background. Figure 5.4 -Comparison of the adaptive approach with [START_REF] Pitié | The linear monge-kantorovitch linear colour mapping for example-based colour transfer[END_REF] (with the post-processing of [START_REF] Rabin | Removing artefacts from color and contrast modifications[END_REF]), [START_REF] Papadakis | A variational model for histogram transfer of color images[END_REF] and [START_REF] Ferradans | Regularized discrete optimal transport[END_REF]. Photo credit: Elektrofisch

Extension to Barycenters

Until now, we focused on the computation of transport between pairs of distributions. Dealing with more that 2 distributions can nevertheless be useful to compute barycenters of color palettes for color normalization of multiple images.

Geodesics between pairs of distributions can indeed be extended to barycenters between an arbitrary finite collection of distributions. Existence and uniqueness of this barycenter is studied in [START_REF] Agueh | Barycenters in the Wasserstein space[END_REF]. Computing the barycenter between discrete distributions requires the resolution of a linear program that corresponds to a multi-marginal optimal transportation, as proved in [START_REF] Agueh | Barycenters in the Wasserstein space[END_REF]. However, in contrast with the case of two distributions, the case of un-weighted sums of Diracs is not anymore equivalent to an assignment problem, which is known to be NP-hard [START_REF] Burkard | Assignment Problems[END_REF]. Computing numerically this barycenter for large scale problems can however be tackled using a non-convex formulation solved with a Lagrangian discretization, which finds applications in image processing [START_REF] Rabin | Wasserstein barycenter and its application to texture mixing[END_REF].

Let us now detail how the previous model can be extended to find a new distribution in-between two or more distributions. Given a set of r = 1

• • • R input histograms µ r ∈ R N associated to clusters Y r ∈ R N ×d ,
and weights (ρ r ) r∈R ∈ (R + ) R , the regularized histogram barycenter µ on the fixed cluster X is a local minimizer of

min µ r∈R ρ r E(µ, µ r ), (5.2)
where

E(µ, µ r ) = min Pr∈Pκ r (µ,µr) κr∈R N ,κr≥0, κr, ν ≥1 C XY , P r + λ G X (D µ P r Y r -X)) 1 + ρ κ r -1 1 . (5.3)
The minimization of (5.2) can be performed by doing a joint minimization on both the barycenter histogram µ and the set of matrices P r and relaxation parameters κ r . This is a non-convex optimization problem. Fortunately, it is separately convex with respect to each of separate group of variables µ and (P r , κ r ) r∈R , so one can use the block coordinate descent algorithm (A. [START_REF] Benamou | A computational fluid mechanics solution to the Monge-Kantorovich mass transfer problem[END_REF]) that iteratively minimizes with respect to each of these variables.

The updates of (P r , κ r ) r∈R can be done in parallel, still using linear solvers. The update of µ is also a linear program which can be solved using for instance interior point solvers [START_REF] Nesterov | Interior Point Polynomial Methods in Convex Programming: Theory and Algorithms[END_REF]. An alternative option is to use first order proximal splitting schemes [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF] as algorithm (A.15), that is well tailored for such highly structured problems. The convergence of such process has been shown in [START_REF] Ferradans | Regularized discrete optimal transport[END_REF].

Remark 2. In the cloud of points formulation where all the mass is uniformly distributed over the clusters X and Y r , the previous problem can be reformulated in terms of clusters positions X instead of histogram weights µ. In the case where λ = 0, ρ → +∞, one recovers the barycenter over the Wasserstein space. As a final example, we show in Figure 5.8 how this method can be applied as a color normalization preprocessing before comparing/registering images of the same object obtained under different illumination conditions.

Experiments

Color normalization is the process of imposing the same color palette on a group of images. This color palette is always somehow related to the color palettes of the original images. For instance, if the goal is to cancel the illumination of a scene (avoid color cast), then the imposed histogram should be the histogram of the same scene illuminated with white light. Of course, in many occasions this information is not available. Following [START_REF] Papadakis | A variational model for histogram transfer of color images[END_REF], an in-between histogram is chosen as the regularized Optimal Transport barycenter and then transferred to each image with the corresponding transport matrix P r .

Finally notice that the entropic regularization of optimal transport can also be considered to compute regularized histogram barycenters in a very efficient way [START_REF] Solomon | Convolutional wasserstein distances: Efficient optimal transportation on geometric domains[END_REF].

Figure 5.8 -The proposed method can be applied as a preprocessing step in a pipeline for objects detection or image registration, where canceling illumination is important. On the first row, a set of pictures of the same object taken at different hours of the day or night is displayed. On the second row, the histogram barycenter has been transfered to the 3 images. Note how the algorithm is able to normalize the illumination conditions of all the images.

Limitations

Due to the one-to-many relaxation and the fact that only the gradient of the average transport flow is penalized, the regularization does not prevent the transport map to associate very different colors to a single pixel or cluster. This leads to undesirable results such as color mixing or color inconsistencies in the modified image. Indeed, new drab colors that do not exist in the target palette ν can be created with regularization. This is illustrated in Figure 5.9, by increasing the regularization parameter λ of model (5.1) or when using instead the entropic regularization of model (4.11).

These experiments show that a large regularization involves an important mixing between colors. Hence, the incorporation of a non-convex constraint that minimizes the variance of colors assigned to each cluster has been proposed in [START_REF] Rabin | Non-convex relaxation of optimal transport for color transfer[END_REF]. ) and the entropic regularization of [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF] are exhibited for color transfer. The second row illustrates the tradeoff between spatial regularity and color fidelity with decreasing regularization parameter λ in (5.1). In the last row, a similar tradeoff is obtained with the method of [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF], when tuning the regularization parameter γ = 1/λ in (4.11).

Non-convex relaxation of regularized Optimal Transport

By using linear programming to optimize the regularized problem (5.1), the dimension of the variables to estimate is greatly increased as one additional variable is needed for each regularization constraint. Simplifications of the regularization term (through the mean transport) is thus needed to reduce the complexity, as illustrated in Figure 4.3.

As previously noticed, such regularization limits the inter-cluster color dispersion but it induces the creation of new drab colors since an important interpolation of the target color palette may occur (i.e. the intra-cluster color variance may be large with the one-to-many assignment). To cope with this issue, the penalization of the dispersion of assigned colors with a non-convex constraint has been introduced in [START_REF] Rabin | Non-convex relaxation of optimal transport for color transfer[END_REF]. With this new formulation, a different optimization tool which decreases the dimension of the problem with respect to linear programming can be considered in order to only deal with the estimation of the transport matrix without considering additional variables for the constraints as in (4.10).

However, to be able to define an algorithm that efficently deals with such non-convex penalization of the color dispersion, differentiable funtional terms have to be considered. Hence, the relaxation here presented is different from the capacity relaxation of the target histogramt proposed in previous section. The closeness to the target histogram is now imposed through a data fidelity term which makes easier the control of the results while simplifying the projection onto the set of acceptable transport matrices.

Dispersion penalization

In order to deal with the aforementioned limitations which are intrinsic to the optimal transport framework, the following relaxed and regularized optimal transport problem can be considered:

P ∈ argmin P ∈Pµ C X,Y , P + ρF ν (P ) + λR(P ) + αD(P ) (5.4)
where:

• P µ = {R N ×N , P ≥ 0, P 1 N = µ} is the convex set of right stochastic matrices (where each row sums to the corresponding bin value in µ),

• F ν (P ) is the fidelity term w.r.t the target histogram ν,

• R(P ) is the regularization prior on the average transport mapping,

• D(P ) measures the dispersion of the average transport map.

To solve efficiently this problem with projected gradient descent, differentiable functions F , R and P are considered.

Linear Constraint The linear constraint P ∈ P µ may be incorporated using an indicator function ι Pµ (P ) = 0 if P ∈ P µ and +∞ otherwise. As the constraint on acceptable transport matrices P µ only concerns the rows, the projection onto this set is simple. The Euclidean orthogonal projector Proj Pµ onto the simplex P µ can be computed in O(nm log(m)) by using, for each row of the matrix, a variant of the Euclidean projector onto the probability simplex (see for instance [START_REF] Chen | Projection onto a simplex[END_REF]).

Fidelity term As the set of acceptable transport matrices P µ does not take into account the target distribution, an additional constraint is required to make sure that the transported histogram P T 1 N is close enough to the target histogram ν. To do so, we rely on the Pearson's χ 2 statistics, which writes for any bistochastic matrix P ∈ P µ

F ν (P ) = 1 2 χ 2 ν (P T 1 n ) = 1 2 1 T n P D 1/2 ν 2 -1 2 , (5.5)
where D µ = diag(ν -1 ) We assume from here, without loss of generality, that ν has non empty bins. Observe that the corresponding fidelity term can be interpreted as a weighted L 2 metric, which further penalizes bins of the target histogram that has small values. This will prevent the model from using very rare features in the exemplar image.

Regularity term A Tikhonov regularization of the gradient of a flow V can be used for incorporating spatial information from the input feature distribution in the gradient operator G X defined on the graph of clusters G X . The norm of the gradient of the mean flow V = D hu P Y -X can therefore be measured as:

R(P ) = 1 2 D µ G X (D µ P Y -X) 2 2 , (5.6)
where the value of the gradient norm is weighted by the corresponding histogram bin value:

(D µ -1 ) ii = µ i .
Dispersion term As pointed out before, the regularized transfer induces a high variability of color assigned to each input color. The simple idea is thus to consider the minimization of the variance of the flow. Denoting Y = D µ P Y , where Y i is the mean color assigned to cluster X i , the intra-cluster variance of the color assigned to X i is defined as Var(Y

) i = (Y 2 ) i -Y 2 i = 1 µ i j P i,j Y j 2 -1 µ i j P i,j Y j 2 .
This variance is finally penalized with respect to each cluster weight to obtain the last functional term:

D(P ) = i µ i Var(Y ) i .
Minimizing this term will promote sparser transport matrices and will prevent strong color mixing in the tackled applications.

Optimization Because of this last term penalizing the variance of the flow, the objective function (5.4) is not convex in P . The inertial Forward-Backward algorithm (A.17) can be used to find a critical point of this non-convex problem [START_REF] Attouch | Convergence of descent methods for semialgebraic and tame problems: proximal algorithms, forward-backward splitting, and regularized gauss-seidel methods[END_REF][START_REF] Ochs | ipiano: Inertial proximal algorithm for nonconvex optimization[END_REF], since it contains a convex non-smooth term (the linear constraint ι Pµ (P )) and the sum of differentiable terms G(P ) = C XY , P + λR(P ) + ρF (P ) + αD(P ), which all have a Lipschitz gradient on the set P µ .

Experiments

Figure 5.10 illustrates the interest of the proposed approach to cancel the effect of transport dispersion. The first row shows input and target histograms, respectively in (a) and (b). The colormap in (a) is arbitrary, and used to illustrate the location reach with the optimal transport map in (b). Observe that the mapping is highly irregular, sending neighboring colors (orange) in two distant locations, as illustrated in (f). Figures (c) and (d), respectively proposed in [START_REF] Ferradans | Regularized discrete optimal transport[END_REF] and [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF], provides more regular mapping of the average flow, as illustrated in (g) and(h), but the colors are spread out to many different locations.

The non-convex model, illustrated in (e), relaxes the matching constraint while controlling the variance of the mapping and the χ 2 statistics with the desired output histogram. This leads to sparser transport flow (i) and almost no color mixing.

The matrices P of transport flows computed with the different methods are illustrated in Figures (f) to (i). The null values P i,j are displayed in white color. These matrices describe how the mass from the bins of the input histogram (in rows) are sent to the bins of the target histogram (in columns). As a bin µ i corresponds to the line i of the matrix, the variance of this bin is larger when a lot of columns j are active (i.e. P i,j > 0) and corresponds to very different features Y j . As illustrated in Figure 5.11, the colors of the target palette are better recovered with increasing penalization of the color variance. When no penalization is applied to the color variance (i.e. α = 0), it corresponds to the model of (5.1). By monitoring the capacity of the target histogram and regularizing the average flow [START_REF] Rabin | Adaptive color transfer with relaxed optimal transport[END_REF], the synthesized images look plausible (5.11b) but they contain new drab colors (that do not exist in the target image) and they are over-smoothed. On the other hand, the transfer is far better when the color variance is penalized (with high values of α). In this case, the final images only contain the colors of the target images.

Discussion

In the future, one interesting aspect to study is the behavior of the non-convex dispersion penalization with respect to the convex entropic regularization. The combination of nonconvex terms with the entropy has for instance been proposed for the domain adaptation problem [START_REF] Courty | Domain adaptation with regularized optimal transport[END_REF]. At first sight, it may seem contradictory to minimize the variance of the flow while maximizing its entropy. Nevertheless, the variance penalization corresponds to the minimization of a local entropy, on the rows of the transport matrix, whereas the method of [START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF] corresponds to the maximization of the global entropy of the matrix. Optimization methods dedicated to problems formulated as the difference of convex functions, known as DC programming methods [START_REF] Horst | Dc programming: Overview[END_REF], may thus be an interesting way to compute efficiently the related transport maps.

Finally, it is important to underline that the regularization of the average optimal transport map allows removing artifacts for color transfer purpose, but its related computational cost is prohibitive and limited in practice to histograms discretized with a few hundreds of bins. On the other hand, the entropic regularization is not designed to deal with transfer artifacts, but it leads to extremely fast computations of approximate transport distances between large scale histograms. Hence, we will see in the next chapter that such regularization is very efficient for involving transport distances in more general problems.

Chapter 6 Application to Image Segmentation

Image segmentation has been the subject of active research for more than 20 years (see e.g. [START_REF] Aubert | Mathematical Problems in Image Processing[END_REF][START_REF] Cremers | A review of statistical approaches to level set segmentation: Integrating color, texture, motion and shape[END_REF] and references therein). For instance, we can refer to the seminal work of [START_REF] Mumford | Optimal approximation by piecewise smooth functions and associated variational problems[END_REF], or to its very popular approximation with level sets developed in [START_REF] Chan | Active contours without edges[END_REF]. This last work provides a very flexible algorithm to segment an image into two homogeneous regions, each one being characterized by its mean gray level value.

In the case of textured images, a lot of extensions of [START_REF] Chan | Active contours without edges[END_REF] have been proposed to enhance the mean value image segmentation model. For instance, local histograms are used in [START_REF] Zhu | Region competition: unifying snakes, region growing, energy/bayes/mdl for multi-band image segmentation[END_REF][START_REF] Ni | Local histogram based segmentation using the wasserstein distance[END_REF], Gabor filters in [START_REF] Vese | A multiphase level set framework for image segmentation using the mumford and shah model[END_REF], wavelet packets in [START_REF] Aujol | Wavelet-based level set evolution for classification of textured images[END_REF] and textures are characterized thanks to the structure tensor in [START_REF] Brox | Unsupervised segmentation incorporating colour, texture, and motion[END_REF][START_REF] Rousson | Active unsupervised texture segmentation on a diffusion based feature space[END_REF].

When considering the global histograms of the regions to segment, there also exist a large range of methods in the literature [START_REF] Aubert | Image segmentation using active contours: Calculus of variations or shape gradients[END_REF][START_REF] Kim | A nonparametric statistical method for image segmentation using information theory and curve evolution[END_REF][START_REF] Brox | Level set segmentation with multiple regions[END_REF][START_REF] Herbulot | Segmentation of vectorial image features using shape gradients and information measures[END_REF]] also based on [START_REF] Chan | Active contours without edges[END_REF]. Recent works make use of the the Wasserstein distance [START_REF] Peyré | Wasserstein active contours[END_REF] to compare globally the histograms. It is important to notice that this class of approaches involves complex shape gradient computations [START_REF] Delfour | Shapes and geometries: analysis, differential calculus, and optimization[END_REF] for the level set evolution equation. Moreover, as these methods all rely on the evolution of a level set function [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on hamilton-jacobi formulations[END_REF], it leads to non-convex methods that are sensitive to the initialization choice and only a local minimizer of the associated energy is computed.

Recently, convexification methods have been proposed to tackle this problem, as in [START_REF] Nikolova | Algorithms for finding global minimizers of image segmentation and denoising models[END_REF][START_REF] Chambolle | An algorithm for mean curvature motion. Interfaces and Free Boundaries[END_REF][START_REF] Chambolle | On total variation minimization and surface evolution using parametric maximum flows[END_REF][START_REF] Berkels | An unconstrained multiphase thresholding approach for image segmentation[END_REF][START_REF] Pock | Global solutions of variational models with convex regularization[END_REF][START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF][START_REF] Brown | Completely convex formulation of the Chan-Vese image segmentation model[END_REF][START_REF] Yıldızoğlu | Active contours without level sets[END_REF]. The original model [START_REF] Chan | Active contours without edges[END_REF] can indeed be convexified, and a global solution can be efficiently computed, for instance with a primal-dual algorithm. Thanks to the coarea formula, a simple thresholding of this global solution provides a global minimizer of the original non-convex problem. Such approaches have not been developed yet for global histogram segmentation with length boundary regularization.

Other models as in [START_REF] Rother | Cosegmentation of image pairs by histogram matching-incorporating a global constraint into MRFs[END_REF][START_REF] Vicente | Joint optimization of segmentation and appearance models[END_REF][START_REF] Ayed | Graph cut segmentation with a global constraint: Recovering region distribution via a bound of the bhattacharyya measure[END_REF][START_REF] Gorelick | Segmentation with non-linear regional constraints via line-search cuts[END_REF][START_REF] Yuan | A fast global optimization-based approach to evolving contours with generic shape prior[END_REF] that use graph-based methods and maxflow formulations [START_REF] Punithakumar | A convex maxflow approach to distribution-based figure-ground separation[END_REF], obtain good minima without level-sets. Nevertheless, these algorithms just consider simple bin-to-bin distances between histograms (for instance 1 or Bhattacharyya) which are not robust enough to deal with coarse quantification of histograms or data outliers.

The use of Optimal Transport for image segmentation has been first investigated in [START_REF] Ni | Local histogram based segmentation using the wasserstein distance[END_REF] for comparing local 1D histograms. In [START_REF] Peyré | Wasserstein active contours[END_REF][START_REF] Mendoza | Linearized multidimensional earth-mover's-distance gradient flows[END_REF], active contours approaches where C µ is the cumulative distribution function of µ and C -1 µ its pesudo-inverse. In the case p = 1, the disance can be simplified as:

W 1 1 (µ, ν) = ||C µ -C ν || 1 .
As noticed in [START_REF] Yıldızoğlu | A convex formulation for global histogram based binary segmentation[END_REF], by defining the M × M cumulative matrix C(i, j) = 1 if i ≤ j and 0 otherwise, the Wasserstein-1 distance between histograms can be included in the model (6.2) by taking:

D(µ, ν) = ||C(µ -ν)|| 1 .
Proximal algorithm such as the Primal-Dual one (A.15) can be thus used to compute a global minimizer u * of the associated non-smooth convex problem. An illustration of the application of the method is given in Figure 6.2. Notice that similar results can be obtained with graph approaches [START_REF] Rother | Cosegmentation of image pairs by histogram matching-incorporating a global constraint into MRFs[END_REF] when 1 norm between 1D histograms are involved. and the M × M matrix C is defined as C i,j = ||F i -F j || 2 , the distance between bins F i and F j .

Approximate MK distance for general histograms

The MK distance can naturally be incorporated into the model (6.2), as proposed in [START_REF] Swoboda | Variational image segmentation and cosegmentation with the wasserstein distance[END_REF]. However, as the MK distance is not smooth, when solving the associated problem with a proximal algorithm, one requires, at each iteration of the algorithm, the computation of the proximal operators of MK(µ 1 1 T u, Hu) and MK(µ 1 1 T (1-u), H(1u)):

Prox MK(µ 1 1 T •,H•) (v) = argmin u 1 2 ||u -v|| 2 + MK(µ 1 1 T u, Hu),
which is a quadratic programming problem that involves inner loops in the process. It has been suggested in [START_REF] Rabin | Convex color image segmentation with optimal transport distances[END_REF] to rely on the entropic regularization of Optimal Transport by extending the works of [START_REF] Cuturi | A smoothed dual approach for variational wasserstein problems[END_REF]. The entropy-regularized Optimal Transport problem is now reformulated as the Sinkhorn distance:

MK λ,≤N (µ, ν) := min

P ∈S(µ,ν) P, C -N λ h(P/N ) + ι H ≤N (µ, ν), (6.5)
where the entropy of the matrix P is defined as h(P ) := -P, log P and the set H ≤N of admissible histograms prescribes the masses of histograms, that are bounded above by N , the number of pixels of the image domain Ω:

H ≤N := µ ∈ R M , ν ∈ R M µ ≥ 0, ν ≥ 0, µ, 1 = ν, 1 ≤ N . (6.6)
By introducing the matrix-valued function Q λ (α, β) i,j = e λ(α i +β j -C i,j )-1 defined for α, β ∈ R M , one can show that the Legendre transform of MK λ,≤N reads:

MK * λ,≤N (α, β) = N λ Q λ (α, β), 1 if Q λ (α, β), 1 ≤ 1 N λ log Q λ (α, β), 1 + N λ if Q λ (α, β), 1 ≥ 1. (6.7)
The constraint ι H ≤N is thus absorbed in this dual formulation. A last interesting point to notice is that this dual function is differentiable and its gradient can be explicitly written as:

∇MK * λ,≤N (α, β) = N (Q λ (α, β)1, Q λ (α, β)1) if Q λ (α, β), 1 ≤ 1 N Q λ (α,β),1 Q λ (α, β)1, Q λ (α, β) T 1 if Q λ (α, β), 1 ≥ 1 . (6.8)
Moreover, it can be shown that the gradient ∇MK * λ,≤N is a Lipschitz continuous function of constant bounded by 2λN .

We emphasis here that we retrieve a similar expression to the one originatively demonstrated in [START_REF] Cuturi | A smoothed dual approach for variational wasserstein problems[END_REF], where the authors consider the Sinkhorn distance on the probability simplex (i.e. the special case where N = 1 and Q, 1 = 1).

Optimization The general final problem we want to solve can be expressed as:

min u MK λ,≤N µ 1 1 T u, Hu + MK λ,≤N µ 0 1 T (1 -u), H(1 -u) + ρ T V (u) + ι [0,1] N (u). (6.9)
Using the Legendre-Fenchel transform, the problem (6.9) can be reformulated as:

min u max α 1 ,β 1 α 0 ,β 0 ,z µ 1 1 T u, α 1 + H u, β 1 + µ 0 1 T (1 -u), α 0 + H(1 -u), β 0 + ∇u, z + ι [0,1] N (u) -MK * λ,≤N (α 1 , β 1 ) -MK * λ,≤N (α 0 , β 0 ) -ι . ≤ρ (z),
where the dual variables α i , β i for the MK term are in R N and the dual variable z for the T V term is in R 2N .

This problem can be efficiently optimized with the Preconditionned formulation [START_REF] Lorenz | An inertial forward-backward algorithm for monotone inclusions[END_REF] of the Primal-Dual algorithm (A. [START_REF] Beck | A fast iterative schrinkage-thresholding algorithm for linear inverse problems[END_REF] proposed in [START_REF] Chambolle | On the ergodic convergence rates of a first-order primal-dual algorithm[END_REF] that combines implicit schemes with the proximal operators of ι [0,1] N (u) and ι . ≤ρ (z) and explicit schemes with the gradient of the other terms. Notice that the process does not require the computation of Optimal Transport and only involves variables of dimensions scaling with N or M . However, the time steps of this algorithm are bounded by the inverse of the Lipschitz constant of the differentiable terms which is proportional to ∇MK * λ,≤N = 2λN . It thus involves a slow explicit gradient ascent in the dual update (A.15) for accurate approximations of the MK distance (i.e. large values of λ). In practice, good results are nevertheless obtained in a few hundred iterations, for relatively high values of λ.

Remark 3. As noticed in [START_REF] Rabin | Convex color image segmentation with optimal transport distances[END_REF], one can resort to the bidualization of MK or MK * λ,≤N in order to remove the dependency of time steps to λ and have algorithms without inner loops. This nevertheless induces a new dual variable which dimension scales with M 2 , the dimension of the transport matrix P .

Experiments

In this experimental section, exemplar regions are defined by the user with scribbles. These regions are only used to built prior histograms, so erroneous labeling is tolerated. Histograms µ 1 and µ 0 are built using hard-assignment on M = 8 d clusters. We use either RGB color features (F = Id and n = d = 3) or the gradient norm of color features (F = ∇. and again n = d = 3). The cost matrix is defined from the Euclidean metric

• in R d space. Final regions are obtained with threshold = 0.5. A few seconds are required to run 100 iterations and segment a 1 Megapixel color image.

Comparison with a local histogram-based approach In order to illustrate the importance of global histogram comparisons over local ones, a local histogram segmentation model is introduced. It looks at the local probability of a pixel to belong to a class (foreground or background). This can be formulated with the functional:

J local (u) = T V (u) + x u(x) h 0 -h V (x) 1 -h 1 -h V (x) 1 , (6.10)
where h V (x) is the histogram estimated on a neighborhood V (x) of x ∈ Ω. Notice that the data term can be computed pointwise (once the map of distances with local histograms has been calculated for the local model). This is a main difference with the previous data term which values depend on the whole state u. On the other hand, it is important to underline that such functionals can be globally minimized on {0; 1} N , using graph methods [START_REF] Boykov | Interactive graph cuts for optimal boundary & region segmentation of objects in N-D images[END_REF] or convexification approaches [START_REF] Chambolle | An algorithm for mean curvature motion. Interfaces and Free Boundaries[END_REF][START_REF] Nikolova | Algorithms for finding global minimizers of image segmentation and denoising models[END_REF][START_REF] Berkels | An unconstrained multiphase thresholding approach for image segmentation[END_REF][START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF].

In Figure 6.3, we can see that the reference color distributions are here not homogeneous in the sense that the histogram of the whole region is not similar to the ones computed in local neighborhoods. Here the orange colors are more probable in the region related to the butterfly. With the local model, the flowers are thus classified as the butterfly, and the darker regions are segmented as being in the background. This Comparison with a shape gradient approach The advantage of the convex model that does not depend on the initialization is finally illustrated through a comparison with the Wasserstein Active Contour method proposed in [START_REF] Peyré | Wasserstein active contours[END_REF].This approach consists in deforming a level set function in order to minimize globally the Wasserstein distance between the reference histograms and the one induced by the segmentation. To make the level set evolve, this formulation requires complex shape gradients computations. As shown in Figure 6.4, even if this model can give good segmentations that are close to the one obtained in Figure 6. General results Figure 6.5 shows the influence of the threshold used to get a binary segmentation. A comparison with the bin-to-bin 1 distance between color histograms is then given in Figure 6.6. This underlines the robustness of Optimal Transport distance with respect to bin-to-bin distance. Contrary to Optimal Transport, when a color is not present in the reference histograms, the 1 distance does not take into account the color distance between bins which can lead to incorrect segmentation.

This robustness is further illustrated in Figure 6.7. It is indeed possible to use a prior histogram from a different image, even with a different clustering of the feature space. The output image u is a regularized weight map that gives the probability of a pixel to belong to the object. This probability map u is finally thresholded with a parameter to segment the input image into a region, which contour is displayed in red. The value = 0.5 is used by default, but other strategies may be defined, such as selecting the threshold value that minimizes the non-relaxed energy (6.1).

Inputs Note that it is not possible with a bin-to-bin metric, which requires the same clustering. Figure 6.8 shows comparisons between the non-regularized model, quite fast but high dimensional model, with the regularized model, using a low dimensional formulation. One can see that setting a large value of λ gives interesting results. On the other hand, using a very small value of λ always yields poor segmentation results.

1 λ = ∞ λ = 1000
Some last examples on texture segmentation are presented in Figure 6.9 where the proposed method is perfectly able to recover the textured areas. We consider in this example the joint histograms of gradient norms on the 3 color channels. Note that the complexity of the algorithm is the same as for color features, as long as we use the same number of clusters to quantize the feature space.

The natural extension of the presented 2 region-based model is the definition of a convex multi-phase segmentation model using Optimal Transport, as done for image restoration in [START_REF] Swoboda | Convex variational image restoration with histogram priors[END_REF]. In this case, the problem of thresholding the final solution to obtain a labeling map may have to be studied more in details. Indeed, even in the case of local

Co-segmentation

As proposed in [START_REF] Swoboda | Variational image segmentation and cosegmentation with the wasserstein distance[END_REF], the previous framework can be extended to deal with the unsupervised co-segmentation of multiple images. The problem of co-segmentation [START_REF] Vicente | Cosegmentation revisited: Models and optimization[END_REF] consists in segmenting simultaneously multiple images that contain the same object of interest without any prior information. In [START_REF] Rabin | Cosegmentation d'images non-supervisée utilisant les distances de sinkhorn[END_REF], the regularized Optimal Transport distance has been considered to solve this problem. Let us consider K images I k defined over domains Ω k , and H k the associated assignment operators that match each pixel of images I k to a bin of the feature space.

The co-segmentation functional defined over the segmentation maps u = (u k ) reads:

J(u) := ι H(u) + K k=1 l>k MK λ H k u k , H l u l + T V (u k ) -α u k 1 (6.11)
where -α u k 1 is a balloon force that prevents from estimating empty areas. In this model, the histograms of regions defined by u k on each image I k are all compared with the MK λ distance. The common mass of all histograms is ensure through the convex set of admissible maps:

H(u) = {u = (u k ), s.t. u k ≥ 0, Ω k u k = constant}.
As illustrated in Figure 6.10, this model does not require any prior and automatically find regions in the different images with common color histograms. One main limitation is that the area of the two regions has to be the same through the constraint ι H , which makes the method not invariant to different scales of the same object. The introduction of unbalanced transport distances [START_REF] Chizat | Unbalanced Optimal Transport: Geometry and Kantorovich Formulation[END_REF] able to compare histograms of different masses is thus the next step to define a more flexible model. Finally notice that the term of functional (6.11) that compares all pairs of histogram with the MK λ distance is linked to the problem of estimating the regularized barycenter between K histograms µ k :

min µ k MK λ (µ k , µ),
which is another interesting aspect to study since regularized barycenters can be computed efficiently with a parallel implementation [START_REF] Solomon | Convolutional wasserstein distances: Efficient optimal transportation on geometric domains[END_REF].

Conclusion of Part II

In this second part, regularized discrete Optimal Transport distances have been presented and applied to color transfer and image segmentation. The aim of the presented works is to enhance the modeling of these image procesing applications using Optimal Transport, while defining fast and unsupervised methods able to automatically process given images.

With the integration of regularization of the transport map, it is possible to design relaxed model able to adapt themselves to heterogeneous color histograms for color transfer purposes. Spatio-color clustering is an important element to decrease the dimension of the Optimal Transport problem, while facilitating the incorporation of spatial consistency in the image domain, and preventing from additional post-processing. The penalization of the color dispersion has been considered with a non-convex model that encourages one-to-one bin associations. Contrary to average transport or entropic regularization of Optimal Transport, it prevents from large mixing of colors contained in the palette of the target image. Notice that the presented methods can be pushed forward to consider other features such that color, gray or gradient patches around each pixel, in order to transfer texture, grain, tone or stylish attributes [START_REF] Efros | Texture synthesis by non-parametric sampling[END_REF][START_REF] Hertzmann | Image analogies[END_REF][START_REF] Pérez | Poisson image editing[END_REF][START_REF] Xia | Synthesizing and mixing stationary gaussian texture models[END_REF][START_REF] Frigo | Motion Driven Tonal Stabilization[END_REF].

It is important to underline that there exist semi-automatic color transfer methods using high level segmentation and/or semantic analysis, while relying on standard Optimal Transport [START_REF] Bonneel | Example-based video color grading[END_REF][START_REF] Frigo | Optimal transportation for example-guided color transfer[END_REF]. They give fast and accurate color transfer results even for videos. However, enhancing the Optimal Transport framework is an important step to make these semi-automatic methods more robust to color outliers.

For the segmentation application, convex methods using user's information and comparing global color distributions have been proposed. The unsupervised automatic cosegmentation of multiple images is a challenge that has been partially addressed. Graph based segmentation approaches can not deal directly with Optimal Transport distances without involving algorithms with inner loops or graph of gigantic size [START_REF] Solomon | Wasserstein propagation for semi-supervised learning[END_REF]. In this context, the dual formulation of the Monge-Kantorivitsch problem is a very relevant solution to propose fast algorithms that directly include the minimization of the transport distances without computing explicitly the transport. Indeed, unavoidable constraints (as the unknown size of the region to segment) can be absorbed within the dual formulation. Moreover, the entropy regularization provides explicit expressions of the gradient of Monge-Kantorovitch distances that are useful for the minimization of general variational problems involving Optimal Transport terms.

General Conclusions and Perspectives

As shown all along the document, generalized Optimal Transport distances including relaxation and regularization can be considered to improve the accuracy performances of image processing applications. New models and algorithms for the continuous and discrete formulations of Optimal Transport have been presented.

Dynamic Optimal Transport

With the fluid dynamic setting, one can compute geodesics between continuous and vanishing densities. The integration of physical regularization of the transport plan is explicit with this formulation and it has made possible the interpolation of ocean images containing complex structures. These works open new possibilities to integrate Optimal Transport distances in numerical weather prediction systems.

Discrete Optimal Transport

In the discrete setting, the regularization of the mean transport plan has been considered for application to color transfer between images. Convex and non-convex models have been proposed in this context in order to define state-of-the-art automatic methods that adapt the proportion of colors contained in the given palette. These methods have also been extended to the computation of barycenters to deal with color normalization of multiple images.

Finally, the entropy regularization of discrete Optimal Transport has been used for image segmentation. A fast and convex model has been designed to segment images, while respecting global color constraints. Whereas it does not reach yet the computational time of global graph-based segmentation approaches, this model is more robust to color histogram quantization thanks to the use of the Optimal Transport distance.

Pros and cons of both formulations With the dynamical formulation, high dimensional continuous and non-negative images can be interpolated and the regularization of the transport map can be done explicitly. Nevertheless, the convergence of this process is slow. On the other hand, the discrete setting can deal with any kind of distributions, but it is limited to small dimensions unless considering the entropic regularization. However, even if such process gives good and fast estimation of the transport cost, it requires positive densities and the transport map is in practice too smooth to be directly used for interpolation purposes.

Future works Problems involving Optimal Transport distances have the potential to impact a very wide range of applications: data assimilation in geophysics [START_REF] Feyeux | Application of optimal transport to data assimilation[END_REF], sound processing, unsupervised clustering through the computation of barycenters [START_REF] Cuturi | Fast computation of wasserstein barycenters[END_REF], statistical analysis of population estimators [START_REF] Bigot | Consistent estimation of a population barycenter in the wasserstein space[END_REF], computer graphics [START_REF] Bonneel | Displacement interpolation using Lagrangian mass transport[END_REF][START_REF] Solomon | Convolutional wasserstein distances: Efficient optimal transportation on geometric domains[END_REF].

Hence, there is a need to provide regularized Optimal Transport distances that can be efficiently minimized using standard optimization algorithms. The recent development of proximal algorithms for non-convex problems (see Appendix A) offers new perspectives for integrating generalized transport distances in a larger range of applications. A specific work is also required to adapt generic models to the targeted applications. Hence, the main perspectives of the presented works will concern the following points.

Geoscience For assimilating satellite images in geosciences, the first direction to investigate is the integration of mass relaxation constraint [START_REF] Piccoli | Generalized wasserstein distance and its application to transport equations with source[END_REF]. There is indeed no reason for an image to have the same mass than the corresponding physical variable in a numerical model. This problem, known as unbalanced Optimal Transport [START_REF] Lombardi | Eulerian models and algorithms for unbalanced optimal transport[END_REF] is crucial for applying the dynamical Optimal Transport framework to real applications, as noticed in [START_REF] Chizat | An Interpolating Distance between Optimal Transport and Fischer-Rao[END_REF] for biological images.

The main limitation of this dynamical approach is the related computational cost. However, when considering adequate basis for decomposing the density and the velocity, faster algorithms can be designed [START_REF] Henry | Optimal transport using Helmholtz-Hodge decomposition and first-order primal-dual algorithms[END_REF]. The integration of the entropic regularization could be a solution to propose accelerated methods but it should include the modeling of advanced physical priors that is necessary to preserve the topology of the structures contained in the images [START_REF] Maas | A generalized model for optimal transport of images including dissipation and density modulation[END_REF].

All these issues are fundamental problems to solve for integrating Optimal Transport in operational data assimilation systems for numerical weather prediction.

Image Processing Other direct perspectives of the works presented in this manuscript are the development of discrete models for transferring general features, multi-phase segmentation models, unbalanced models [START_REF] Frogner | Learning with a wasserstein loss[END_REF][START_REF] Chizat | Unbalanced Optimal Transport: Geometry and Kantorovich Formulation[END_REF] for co-segmentation or the integration of additional priors such as the shape of the objects [START_REF] Schmitzer | Modelling convex shape priors and matching based on the gromov-wasserstein distance[END_REF][START_REF] Schmitzer | Globally optimal joint image segmentation and shape matching based on wasserstein modes[END_REF].

The interest of Optimal Transport is nevertheless not limited to image segmentation or color transfer and its relation with image denoising or cartoon decomposition for instance has already been shown [START_REF] Lellmann | Imaging with kantorovich-rubinstein discrepancy[END_REF]. A yet not studied problem concerns its use in non-local models. The introduction of Optimal Transport for defining regular non-local weight distributions between adjacent pixels is for instance an option to enhance the construction of nonlocal graphs.

Computer Vision and Machine Learning Various works on color image retrieval [START_REF] Pele | A linear time histogram metric for improved sift matching[END_REF] or artistic image indexation [START_REF] Hurtut | Adaptive image retrieval based on the spatial organization of colors[END_REF] have shown on small datasets that the Optimal Transport framework is able to increase retrieval performances. Spatial information is typically only considered to enhance the feature description [START_REF] Sharma | Expanded parts model for human attribute and action recognition in still images[END_REF] but not used in the metric itself. Partial transport models, that relax the constraint of transporting the whole histograms mass, are nevertheless still required to increase the robustness to data outliers.

The development of approximate Monge-Kantorovitch distances with the entropic regularization has made possible the real-time comparison of histograms quantized with thousands of bins. In order to transfer such tools for data retrieval, a research effort is still required. It involves the learning of the ground cost [START_REF] Cuturi | Ground metric learning[END_REF] inside the transport distance as well as efficient implementation, to reach the performances of local descriptors and more simple metrics in optimized frameworks [START_REF] Jégou | Aggregating local image descriptors into compact codes[END_REF].

Data analysis and Clustering

The definition and computation of barycenters relative to Optimal Transport distance is a major issue for proposing clustering algorithms adapted to the analysis of structured data. The notion of mean and variance with respect to the Optimal transport distance is then necessary. Recent studies on Wasserstein barycenters [START_REF] Agueh | Barycenters in the Wasserstein space[END_REF] have led to practical algorithms both for standard [START_REF] Carlier | Numerical methods for matching for teams and Wasserstein barycenters[END_REF] or regularized barycenters [START_REF] Ferradans | Regularized discrete optimal transport[END_REF] [START_REF] Solomon | Wasserstein propagation for semi-supervised learning[END_REF].

An interesting direction of research is the definition of geodesic Principal Component Analysis in the Wasserstein space [START_REF] Bigot | Geodesic PCA in the Wasserstein space[END_REF]. First algorithmic works [START_REF] Seguy | An Algorithmic Approach to Compute Principal Geodesics in the Wasserstein Space[END_REF] have been developed with applications to color transfer, but this subject deserves to be explored more in depth to propose new clustering algorithms.

Audio Processing The transfer of Optimal Transport tools into the audio community is an open subject of research. As for image processing, the transport distance could define a robust dissimilarity measure and the transport map could be used for interpolation purposes (voice, prosody, music, etc).

A fundamental difference with images is that the audio signal is described with phase and modulus, and can not be represented by non negative densities to apply the Optimal Transport framework. Computing transport distances between distributions of features extracted from audio signal is nevertheless possible. Dictionary learning models could also be considered [START_REF] Cuturi | A smoothed dual approach for variational wasserstein problems[END_REF] to separate different sources presents in a signal. For audio interpolation, the use of the Optimal Transport map is more challenging and dedicated tools should be defined. Works on phase retrieval from signal modulus [START_REF] Mallat | Phase retrieval for the Cauchy wavelet transform[END_REF] are thus an interesting direction to handle such issue.

In conclusion This document presents contributions in the field of Optimal Transport for image processing that have been realized thanks to the help of a lot of colleagues that I would like to thank: Jean-François Aujol, Afaf Bouharguane, Vincent Chabot, Sira Ferradans, Romain Hug, Emmanuel Maitre, Quentin Mérigot, Maëlle Nodet, Edouard Oudet, Gabriel Peyré, Julien Rabin, Arthur Vidard and Romain Yıldızoğlu. I acknowledge Antonin Chambolle, Patrick Pérez and Christoph Schnörr for accepting to review this thesis as well as Jean-François Aujol, Jérémie Bigot, Julie Delon, Angelo Iollo and Étienne Mémin for agreeing to be committee members.

In the following, we will consider smooth and non-smooth functions. We will refer to a smooth function as soon as it can be differentiated once ∀x ∈ E, i.e. f ∈ C k (E, R), for k ≥ 1. A smooth functional is thus differentiable and its gradient will be denoted as ∇f . Hence, a function f is called non-smooth if is at most in C 0 (E, R). Definition 6 (Lipschitz). A function f defined from E to R ∪ +∞, is L-Lipschitz is there exists L > 0 such that ∀(x, y) ∈ E × E:

||f (x) -f (y)|| ≤ L||x -y||.
if L < 1 then f is said firmly non-expansive.

In order to minimize a convex and differentiable function that admits minimizers, gradient descent algorithms are standard optimisation tools. Second-order methods such as the Newton's one, can be considered to algorithms with faster convergence rates. Before detailing the non-smooth context, we finally recall the explicit first order gradient descent algorithm.

Theorem 2 (Gradient descent Algorithm). Let f be a convex, proper, l.s.c and differentiable function defined from E to R ∪ +∞ with a L-Lipschitz gradient. Let τ be a real parameter checking τ < 1 L , so that the operator Id -τ ∇f is firmly non-expansive. Let

x 0 ∈ E, then the sequence (x n ) n∈N defined ∀n ∈ N as

x n+1 = (Id -τ ∇f )(x n ), (A.1)
converges to a minimizer of f .

A.1.1 Subdifferential

In order to consider the minimization of non-smooth functionals, we now recall standard results.

Definition 7 (Subdifferential). Let f be a convex function defined from E to R ∪ +∞, the subdifferential of f at point x ∈ E is the multivalued operator that associates to x the set: ∂f (x) = {u ∈ E such that ∀y ∈ E, y -x, u + f (x) f (y)}.

The subdifferential of a function at a point x can be empty or it can be a convex set. If f is convex et differentiable ∀x ∈ E, the subdifferential is reduced to the singleton {∇f (x)}, the gradient of f at point x.

The minimizers of non-smooth convex functions can be characterized with the subdifferential.

Theorem 3. Let f be a convex, l.s.c and proper function defined from E to R ∪ +∞, then argmin f = zer ∂f = {x ∈ E such that 0 ∈ ∂f (x)}.

The subdifferential also has the summability property:

Lemma 1. Let f and g be two convex, l.s.c and proper functions. If dom (f ) ∩ dom (g) = ∅ then, for J = f + g, we have ∂J = ∂f + ∂g.

A.1.2 Proximal operator

The proximal operator associated to a convex function is defined as follows.

Definition 8 (Proximal operator). Let f be a convex function and x ∈ E, Prox f (x) is the unique minimizer of the strictly convex function:

x → f (y) + 1 2

x -y 2 .

(A.

2)

The proximal operator can be seen as the generalization of the concept of the projection onto a convex set. Indeed, if f is the characteristic function of a close convex set C then Prox f (x) is the orthogonal projection of x onto C.

In the following, we will say that f is simple if its proximal operator can be computed explicitly.

When f is not a convex function; then the proximal operator may not be defined or it can be multivalued. The convexity of f is a sufficient but not necessary condition for the existence and uniqueness of Prox f (x). In the following, a function will be said proximable if the proximal operator Prox f (x) is defined for every point x ∈ E. A proximable function can therefore be non-convex.

An important property of the proximal operator is now given. Proposition 3. Let f be a convex, l.s.c and proper function defined on E, then Prox f is a 1-Lipschitz operator.

With the proximal operator, it is thus possible to minimize non-smooth convex functional through an implicit descent method.

Theorem 4 (Proximal Algorithm). Let f be a convex, l.s.c and proper function defined from E to R ∪ +∞, then the fixed points of Prox f are minimizers of f . Hence, ∀τ > 0 and x 0 ∈ E, the sequence (x n ) ∈N defined as

x n+1 = Prox τ f (x n )
converges to a minimizer of f .

A.1.3 Conjugate function

A large range of approaches use the duality of convex functionals to design minimization algorithms. For all function f from E to R, one can associate a conjugate function f * through the following relation. (A.3)

Conjugate functions plays an essential role in the optimisation of convex functionals. There exists strong relations between conjugation subdifferential and proximal operators, from which the equivalence between primal and dual problems can be established. One main property of this transform is the following.

Theorem 5 (Biconjugate). All function f defined on E with values in [-∞, +∞], convex, l.s.c and proper is equal o its biconjugate: f * * = f . One useful corollary can be shown from this relation.

Corollary 1 (Moreau's identity). If f is convex, l.s.c and proper, then , for all τ > 0: Prox τ f (x) + τ Prox f * /τ (x/τ ) = x, that is to say Prox f + Prox f * = Id.

A.2 Convex functionals

Standard algorithms for the minimization of convex and non-smooth functionals are now introduced. We will assume that the considered functions admit at leat one minimizer.

A.2.1 Primal Proximal Algorithms

In this subsection, we first consider the problem:

min x∈E J(x) = f (x) + g(x), (A.4)
where f and g are convex, l.s.c, proper, coercive and bounded from below functions.

We first assume that f is differentiable with a L-Lipschitz gradient and that the condition of Lemma 1 is satisfied so that the subdifferential of J is equal to the sum of the subdifferentials of f and g for all point of the domain of ∂J. From Theorem 3, a minimizer x * of J then checks -∇f (x * ) ∈ ∂g(x * ). Under these assumptions, a minimizer of J can be estimated with the Forward-Backward (FB) Algorithm.

Theorem 6 (Forward-Backward Algorithm). Let J = f + g be the sum of two convex, l.s.c, proper and coercive functions. Let f be differentiable with a L-Lipschitz gradient.

Let τ < 1 L and x 0 ∈ E then the sequence (x n ) n∈N defined ∀n ∈ N as

x n+1 = Prox τ g (Id -τ ∇f )(x n ) (A.5)
converges to a minimizer of J.

A complete description of this algorithm can be found in [START_REF] Combettes | Proximal splitting methods in signal processing[END_REF]. Notice that using an sadequate sequence of time steps τ n , it is possible to accelerate its convergence [START_REF] Beck | A fast iterative schrinkage-thresholding algorithm for linear inverse problems[END_REF][START_REF] Chambolle | On the convergence of the iterates of the "fast iterative shrinkage/thresholding algorithm[END_REF].

where K is defined on E with values in F , f defined from F to R and g defined from E to R. We now assume that the proximal operator of f • K can not be computed explicitly, whereas f and g are simple (i.e. their proximal operators can be computed explicitly). In order to define a simple algorithm for minimizing J in this case, one can consider DR or ADMM algorithms by introducing an auxiliary variable y ∈ F and solve the following problem: By denoting f = f + g and g = ι K , the Douglas-Rachford Algorithm can indeed be directly applied to the function f + g with the variable z = (x, y). The proximal operator of f (z) can be computed separately for the two variables x and y through the proximal operators of f and g. The proximal operator of g(z) can be computed explicitly as the orthogonal projection onto the linear constraint Kx = y.

The Primal-Dual formulation is another way to deal with such problem. The convex conjugate (A.3) can indeed be used to separate the function f from the operator K: where K * denotes the adjoint operator of K defined from F to E. Injecting this expression in the problem (A.10), ones gets a new problem where the function f has been dualized. As soon as the subdifferential of f * is no empty for all point of its domain, the above supremum is a maximum. This will be the case for all the functionals we will consider, so that we can turn into the general saddle-point problem (x * , y * ) ∈ E × F : Notice that problems (A.10) and (A.13) are equivalent in the sense that if x is a solution of (A.10) then there exists y such that the pair (x, y) is solution of (A.13). Conversely, if (x, y) is a solution of (A.13) then x is a solution of (A.10).

In [START_REF] Chambolle | A first-order primal-dual algorithm for convex problems with applications to imaging[END_REF], a variant of the Arrow-Hurwicz Algorithm [START_REF] Arrow | Studies in linear and non-linear programming[END_REF] has been proposed to deal with this problem. Such algorithm is equivalent to the DR one, if K = Id. We now recall a more general expression of this algorithm [START_REF] He | Convergence analysis of primal-dual algorithms for a saddlepoint problem: from contraction perspective[END_REF][START_REF] Vũ | A splitting algorithm for dual monotone inclusions involving cocoercive operators[END_REF][START_REF] Condat | A primal-dual splitting method for convex optimization involving lipschitzian, proximable and linear composite terms[END_REF][START_REF] Bot | On the convergence rate improvement of a primal-dual splitting algorithm for solving monotone inclusion problems[END_REF][START_REF] Chambolle | On the ergodic convergence rates of a first-order primal-dual algorithm[END_REF][START_REF] Lorenz | An inertial forward-backward algorithm for monotone inclusions[END_REF][START_REF] Drori | A simple algorithm for a class of nonsmooth convex-concave saddle-point problems[END_REF], able to deal with an additional differentiable convex function h: (A.14)

Theorem 9 (Primal-Dual Algorithm). Let f , g and h be three convex, l.s.c and proper functions, f defined from F to [-∞, +∞], g defined from E tp [-∞, +∞], h differentiable defined from E to [-∞, +∞] with a L-Lipschitz gradient and K a linear operator defined from F to E. Let L K = K , τ > 0 and σ > 0 such that (1/τ + L)/σ > L 2 K .

Theorem 10 (Inertial Forward-Backward Algorithm). Let f and g be two l.s.c and proper functions. Let f be differentiable with a L-Lipschitz gradient, g be proximable and J = f + g verifying the KL property. Let β ∈ [0; 1[, τ < 2(1 -β)/L and x 0 ∈ E, then the sequence (x n ) n∈N defined as

x n+1 = Prox τ g (x n -τ ∇f (x n )) + β(x n -x n-1 ) (A.17)

converges to a critical point x of J.

The Forward-Backward can also be extended to deal with several proximable terms i g(x) [START_REF] Chouzenoux | A block coordinate variable metric forward-backward algorithm[END_REF]. Finally notice that if f is non-convex but c-semi convex (such that f (x) + c||x|| 2 is convex), other proximal algorithms in the primal [START_REF] Kowalski | Thresholding rules and iterative shrinkage/thresholding algorithm: A convergence study[END_REF] or primal-dual [START_REF] Möllenhoff | The primal-dual hybrid gradient method for semiconvex splittings[END_REF] formulations also ensure the convergence to a critical point of J.

A.3.2 Non-convex coupling

For some particular applications, we finally consider the following problem: where f and g are two convex functions and h is non-convex but differentiable. We make the additional assumption that J is separately convex in x and y but not in (x, y).

As the function h may not be semiconvex, the abovementioned approaches can not be applied. The estimation of a critical point of such functional is ensured by an alternate minimization by blocks [START_REF] Tseng | Convergence of a block coordinate descent method for nondifferentiable minimization[END_REF]:

Theorem 11 (Block Coordinates Descent Algorithm). Let J(x, y) = f (x) + g(y) + h(x, y) such that for all fixed (x, ỹ) ∈ E × E , J(., ỹ) et J(x, .) are convex, l.s.c and proper functions. Let (x 0 , y 0 ) ∈ E × E , then the sequence (x n , y n ) n∈N defined by converge to a critical point (x * , y * ) of J.

This approach nevertheless requires to solve sub-optimization problems at each iteration of the algorithm. For some very specific functionals (A.18), these intern loops can be avoided and primal-dual approaches can be considered [START_REF] Pierre | Luminancechrominance model for image colorization[END_REF].

In fact, as soon as f and g are proximable functions (even if they are not convex), h is a C 2 function and the gradient ∇ x h(., y) for y fixed is L y -Lipschitz (resp. ∇ y h(x, .) is L x -Lipschitz), a proximal alternating algorithm can be applied [START_REF] Bolte | Proximal alternating linearized minimization for nonconvex and nonsmooth problems[END_REF]. This last approach is related to the application of the Forward-Backward algorithm on the problem f (z) = f (x) + g(y) and g(z) = h(z) with the variable z = (x, y).

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Limitations of Optimal Transport in Imaging . . . . . . . . . . . . . . . . 1.3 Contributions and Organization of the document . . . . . . . . . . . . . I Continuous Optimal Transport 2 Dynamic Optimal Transport 2.1 Estimation of the Transport Map . . . . . . . . . . . . . . . . . . . . . . 2.2 Fluid Mechanics Formulation . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Convexification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Discretization and Optimization . . . . . . . . . . . . . . . . . . . . . . 2.5 Limitations and Motivations . . . . . . . . . . . . . . . . . . . . . . . . . 3 Generalizations for Image Interpolation 3.1 Generalized Cost Functions . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.1 Interpolation between L 2 -Wasserstein and H -1 . . . . . . . . . . . 3.1.2 Modelling of obstacle . . . . . . . . . . . . . . . . . . . . . . . . . 3.1.3 Anisotropic transport . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Optimal Transport with physical priors . . . . . . . . . . . . . . . . . . . 3.2.1 Non-convex coupling . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.2 Synthetic tests . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.3 Image interpolation in oceanography . . . . . . . . . . . . . . . .

Figure 1 . 1 -

 11 Example of massive use of images. (a) Reconstruction of the Colosseum from Flickr images [2]. (b) UrbanDive navigation system of IGN.

Figure 1 . 2 -

 12 Figure 1.2 -Different examples of satellite images used for data assimilation. (a) Altimetric reconstruction from JASON satellite data. (b) Ocean color/ Chlorophyll from the MODIS captor of ENVISAT satellite. (c) Sea surface Temperature from the MODIS captor of ENVISAT satellite.

Figure 1 . 4 -

 14 Figure 1.4 -Illustration of the limitation of Optimal Transport for color transfer. The Optimal Transport map T between the color distributions of the middle parrot µ 0 and the left parrot µ 1 is computed. The map is then used to modify the colors of the middle parrot in order to obtain the right one. As the transport map is highly irregular, very different colors are assigned to pixels of the background that were originally close both in spatial and color spaces.

Figure 1 . 5 -

 15 Figure 1.5 -Illustration of the transport map in the (a) continuous and (b) discrete formulations of the Optimal Transport problem. In (b), white color represents P i,j = 0.

Figure 2 . 1 -

 21 Figure 2.1 -Transport between characteristic functions. Evolution of ρ (•, t) for several values of t. The red curve denotes the boundary of the area with positive density.

Figure 2 . 2 -

 22 Figure 2.2 -Illustration of the convergence of the minimization process. (a) The cost function J reaches its minimum with a thousand iterations. (b) The convergence of the iterates ρ k to ρ * is slower. The reference solution ρ * has been previously obtained with 10 8 iterations. Different proximal algorithms (ADMM, symmetric and asymmetric Douglas Rachford, Primal Dual) tunned with their best parameters are compared.

Figure 2 . 4 -

 24 Figure 2.4 -Two Gaussians examples. The Gaussians, represnted by their level lines, split along the optimal path between ρ 0 and ρ 1 .

Figure 3 . 1 ,

 31 which presents the level-lines of the estimated densities ρ(t, •) for A = Id d and different values of β. It shows the evolution of the solution between a linear interpolation of the densities (β = 0) and a displacement interpolation with transport (β = 1).

Figure 3 . 1 -

 31 Figure 3.1 -Display of the level sets of ρ(t, •) for several value of t and β. For t = 0 and t = 1, ρ(t, •) exactly corresponds to ρ 0 and ρ 1 .

Figure 3 . 2 -

 32 Figure 3.2 -Interpolation between Gaspard Monge (ρ 0 ) and Leonid Kantorovitch (ρ 1 ). The first line is the same as in Figure 2.3. On the second line, by taking β = 0.75, a smoother interpolation between images is estimated.

Figure 3 . 3 -

 33 Figure 3.3 -Evolution of ρ(t, •) for several values of t, using a Riemannian manifold with weights w(x) (constant in time) restricting the densities to lie within a 2-D static labyrinth map.

  ure 3.4), there are no equivalence with any static Optimal Transport problem and the proof of existence of solutions remains an open subject.

Figure 3 . 4 -

 34 Figure 3.4 -Evolution of ρ(t, •) for several values of t, using a Riemannian manifold with weights w(t, x) (evolving in time) restricting the densities to lie within a 2-D dynamic labyrinth map (i.e. with moving walls in green and red).

Figure 3 . 5 -

 35 Figure 3.5 -Illustration of the mass ρ * (t) estimated between ρ 0 and ρ 1 , through the computation of the transport costs defined by two different anisotropic domains A k , illustrated by the blue directions in the first column of the two last rows. The first row is the isotropic transport.

= ρ 1 Figure 3 . 6 -Figure 3 . 7 -

 13637 Figure 3.6 -Two Gaussians experiments with penalization and a null initialization.Plot of the isolevels of the density ρ(t) along the different computed optimal path. The top line is realized with incompressibility (i.e. divergence free) penalization, the second with a translation penalization and the last one with a rigid penalization.

= ρ 1 Figure 3 . 8 -Figure 3 . 9 -

 13839 Figure 3.8 -Bar experiment.Plot of the isolevels of the density ρ(t) along the optimal path computed with different approaches. The first line is the Optimal Transport, the second is the proposed approach with an incompressible penalization and the last one with a rigid penalization that conserves the nature of the object to transport.

1 Figure 3 . 10 -

 1310 Figure 3.10 -Interpolation of Sea Surface Height images (ρ 0 and ρ 1 ) in oceanography: the creation of vortexes in Cap Point is finely estimated.

Figure 4

 4 Figure 4.1 -(a) Optimal Transport matrix P between histograms X (in blue) and Y (in red). The white color correspond to P i,j = 0. (b) The transport map does not impose one-to-one mappings betweens bins and the bin X i is transported to several bins Y j .

  κ = (1, 1, 0.1, 1.5) κ = (0, 2, 1, 1) κ = (0.1, 10, 0.1, 10)

Figure 4 . 2 -

 42 Figure 4.2 -Relaxed transport computed between X (blue dots) and Y (red dots)for different values of κ. Note that κ = (1, 1, 1, 1) corresponds to classical Optimal Transport. The mappings P i,j that relate X i and Y j are plotted as line segments. The segments are dashed if P i,j ∈]0.1, 1[ and solid if P i,j = 1. The parameter κ = (1, 1, 0, 2) mean that the sum of links connecting each red point to blue points equals 1, while the sum of links connecting each blue point to red points is in [0; 2]. In this case, as all the parameters are integers, assignments are obtained (solid lines) and each red point is assigned to exactly one blue point, whereas each blue point can be connected to 0, 1 or 2 red points.

Figure 4 . 3 -

 43 Figure 4.3 -Computation of the mean transport map T from the Optimal Transport matrix P between histograms X (in blue) and Y (in red). (a) The bin X i is transported to several bins Y j . (b) The mean transport from bin X i is T (X i ) = Ȳi , obtained from relation (4.8). The white color correspond to P i,j = 0.

Figure 4 . 4 -

 44 Figure 4.4 -Given two sets of points X (in blue) and Y (in red), we show the points Z = diag(P 1) -1 P Y (in green), and the mappings P i,j as line segments connecting X i and Y j , which are dashed if P i,j ∈]0.1, 1[ and solid if P i,j = 1. The results were obtained with the relaxed and regularized Optimal Transport formulation, setting the parameters to κ = (0.1, 8, 0.1, 8). Note the influence of a change in λ X and λ Y on the final result: with no regularization (λ X = λ Y = 0) only few points in the data set are matched. The introduction of regularization (λ X = λ Y = 0.001) spreads the connections among the clusters, while maintaining the cluster-to-cluster matching. For a high value of λ X = λ Y = 10, the regularization tends to match the clusters with similar shape with each other, where the shape is defined by the graph structure. The graphs G X and G Y are represented with the nodes on blue and red respectively, and the edges as solid lines.

Figure 4 . 5 -

 45 Figure 4.5 -Illustration of the effect of different regularizations of Optimal Transport computed between histograms (a) and (b). The standard Optimal Transport map is shown in (c). Relaxaed and Regularized Optimal Transport flow of[START_REF] Ferradans | Regularized discrete optimal transport[END_REF] is presented in (d) and the entropic reguarization of[START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF] in (e).

Figure 5 . 1 -

 51 Figure 5.1 -Artifacts of color transfer with Optimal Transport. (a) Input image I x . (b) Target color palette I Y . (c) Exact color specification with optimal transport. (d) Post-processing with spatial regularization.

Figure 5 .Figure 5 . 2 -

 552 Figure 5.2 -Example of color transfer. Given images I X and I Y with their corresponding 3-D RGB color distributions µ X and µ Y (represented here using the 2-D projection of every pixel in the RG plane), the goal is to define an image I Z that has the geometry of I X and a histogram µ Z that is similar to µ Y .

  ), and similarly for the measure ν obtained from image I Y .

Figure 5 . 3 ,Figure 5 . 3 -

 5353 Figure 5.3 -Color transfer process. Super-pixel segmentation of (a) source image I X and (b) target image I y . Color transfer without (c) and with (d) relaxation (i.e. ρ = ∞ and ρ = 10 3 ). Photo credit: Nicolas Le Dilhuit

Figure 5 .

 5 6 finally presents examples of color transfer on HDR images.

Figure 5 . 5 -Figure 5 . 6 -

 5556 Figure 5.5 -Adaptive color transfer with relaxed Optimal Transport. (a) Input image I x . (b) Target color palette I Y . (c) Transfer with adaptive relaxation. (d) Relaxed color palette κ ⊗ ν (the darker, the less a color is used). The model adapts itself to the color distributions and automatically selects the adequate proportion of color to transfer.

Figure 5 .Figure 5 . 7 -

 557 Figure 5.7 -Results for the barycenter algorithm on different images. Note how as ρ approaches (0, 1), the histogram of the barycenter image becomes similar to the histogram of I Y .

Figure 5 . 9 -

 59 Figure5.9 -Creation of drab colors with the regularization. The influence of the regularization of the average transport (5.1) and the entropic regularization of[START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF] are exhibited for color transfer. The second row illustrates the tradeoff between spatial regularity and color fidelity with decreasing regularization parameter λ in (5.1). In the last row, a similar tradeoff is obtained with the method of[START_REF] Cuturi | Sinkhorn distances: Lightspeed computation of optimal transport[END_REF], when tuning the regularization parameter γ = 1/λ in (4.11).

Figure 5 . 10 -

 510 Figure 5.10 -Illustration with 1D histograms of various models for regularized optimal transport.

  (a) Input (b) α = 0 (c) α = 10 (d) α = 100 (e) Exemplar

Figure 5 . 11 -

 511 Figure 5.11 -Color transfer. The colors of the exemplar images (e) are transfered to the input images (a).

Figure 6 . 2 -

 62 Figure 6.2 -Segmentation of a natural image. The input histogram prior µ 1 (resp. µ 0 ) is given by the white (resp. green) area. The final contour is obtained by thresholding the solution with = 0.5.

For, P 1 M

 1 general feature space F of dimension d > 1, the previous reformulation in term of cumulative histograms does not hold. Hence, we would like to consider D as the Monge-Kantorovitch distance between the histograms µ, ν ∈ R M : MK(µ, ν) := min P ∈S(µ,ν) of admissible transport matrices is S(µ, ν) := {P ∈ R M ×M + = µ and P T 1 M = ν}. (6.4)

Figure 6 . 3 -

 63 Figure 6.3 -Segmentation of an image where local histograms (9 × 9 neighborhoods) are different from the global ones. The user defines scribbles which indicate the object to be segmented (here in white) and the background to be discarded (in green). The segmentation fails for the local histogram model as it classifies the orange areas in the first class and the darker ones in the second class.

Figure 6 . 4 -

 64 Figure 6.4 -The Wasserstein active contours method [165] has been initialized in two different ways (a-b), the corresponding segmentations being presented in (c-d). When carefully initialized, it leads to a segmentation close to the one obtained with the convex approach (see Figure 6.3 (c)).

Figure 6 . 5 -

 65 Figure6.5 -Illustration of the influence of the threshold parameter on the segmentation result. The output image u is a regularized weight map that gives the probability of a pixel to belong to the object. This probability map u is finally thresholded with a parameter to segment the input image into a region, which contour is displayed in red. The value = 0.5 is used by default, but other strategies may be defined, such as selecting the threshold value that minimizes the non-relaxed energy (6.1).

Figure 6 . 6 -

 66 Figure 6.6 -Comparison of the segmentation results obtained from the segmentation models using MK λ distances, together with the bin-to-bin 1 distance. The same regularization parameter ρ is used for every segmentations. Note that the Optimal Transport similarity measure is a more robust statistical metric between histograms than 1 .

Figure 6 . 7 - 10 Figure 6 . 8 -

 671068 Figure 6.7 -Illustration of the interest of Optimal Transport for the comparison of histograms. Its robustness makes it possible to use prior histograms from different images (here histograms are taken from image 1 and used to segment images 2 and 3), even with a different clustering of feature space. Note that it is not possible with bin-to-bin metric, which requires the same clustering.

Input λ = ∞ Figure 6 . 9 -

 ∞69 Figure 6.9 -Texture segmentation using joint histograms of color gradient norms.

u 1 u 2 Figure 6 . 10 -

 2610 Figure 6.10 -Example of co-segmentation of two images.

Definition 9 (

 9 Convex conjugate). Let f be a function defined from E to R, the Fenchel conjugate of f is the convex function defined from E to R byf * (u) = sup x∈E ( x, u -f (x)).

  min x∈E f (Kx) + g(x) = min x∈E,y∈F f (y) + g(x) + ι K (x, y), (A.11)where h is the characteristic function of the linear constraint Kx = y:ι K (x, y) = 0 if Kx = y +∞ otherwise.

  ∀x ∈ E, f (Kx) = sup y∈F ( y, Kx -f * (y)) = sup y∈F ( K * y, x -f * (y)) (A.12)

K

  * y, x -f * (y) + g(x). (A.[START_REF] Ayed | Graph cut segmentation with a global constraint: Recovering region distribution via a bound of the bhattacharyya measure[END_REF]

K

  * y, x -f * (y) + g(x) + h(x).

  min x∈E,y∈E J(x, y) = f (x) + g(y) + h(x, y), (A.18)

f

  (x) + h(x, y n ), y n+1 = argmin y g(y) + h(x n+1 , y), (A.19)

In[START_REF] Maas | A generalized model for optimal transport of images including dissipation and density modulation[END_REF], a different choice has been made and the problem is solved on (ρ, v) while penalizing the non convex constraint ||∂ρ + div(ρv)||

.

Part I Continuous Optimal Transport

Convex histogram-based image segmentation using the Wasserstein distance for comparing global multi-dimensional histograms of the region of interest have then been proposed. Again, such non-convex methods are sensitive to the initial contour. Moreover, their computational cost is too important (several minutes) even if they include some approximations of the Wasserstein distance.

In order to include more robust metrics between histograms while being independent of any initialization choice, we now describe some convex [START_REF] Swoboda | Variational image segmentation and cosegmentation with the wasserstein distance[END_REF] and fast [START_REF] Yıldızoğlu | A convex formulation for global histogram based binary segmentation[END_REF][START_REF] Rabin | Convex color image segmentation with optimal transport distances[END_REF][START_REF] Rabin | Cosegmentation d'images non-supervisée utilisant les distances de sinkhorn[END_REF] segmentation models including Optimal Transport between global color histograms.

Convex histogram-based image segmentation

Let I : x ∈ Ω → I(x) ∈ R n be a color image, defined over the N -pixels domain Ω (N = |Ω|), and F a feature-transform of n-dimensional descriptors (i.e. colors, gradient orientations...) such that FI(x) ∈ R d . The feature space F is assumed to be clustered in M regular bins F j ∈ R d . Hence, one can design a linear assignment operator H of size M × N such that for each pixel x, H(x, j) = 1 if FI(x) ∈ F j and 0 otherwise.

The objective is to segment the image I into 2 regions, by defining a binary segmentation map u : Ω → {0, 1} of the whole image domain with respect to two fixed probability distributions of features µ 1 for the foreground and µ 0 for the background defined in F . This prior information is commonly obtained by scribbles given by a user in an interactive segmentation process [START_REF] Boykov | Interactive graph cuts for optimal boundary & region segmentation of objects in N-D images[END_REF]. In a local model, the distance to the scribble information is considered to perform the segmentation, using for instance local histogram comparison [START_REF] Blake | Interactive image segmentation using an adaptive gmmrf model[END_REF][START_REF] Ni | Local histogram based segmentation using the wasserstein distance[END_REF] or geodesic distances [START_REF] Criminisi | Geodesic image and video editing[END_REF]. As illustrated in Figure 6.1, the goal is here different as global statistics of the regions to segment are compared. As proposed in [START_REF] Yıldızoğlu | A convex formulation for global histogram based binary segmentation[END_REF], the following model is considered:

where T V (u) = |Du| is the total variation of u measuring the length of its level sets and ρ ≥ 0 is the regularization parameter. The function D is a convex distance between feature distributions. Finally, π(u) is the empirical discrete probability distribution of features FI using the binary map u which writes:

where |u| 1 = x |u(x)| = u, 1 = 1 T u for positive values of u denotes the size of the region defined by u = 1, i.e. the mass of the histogram Hu. The functional (6.1) is highly non-convex since π is a non linear operator, and the objective is to find a minimum over the non-convex set {0, 1} N .

In order to avoid the use of combinatorial algorithms for solving this problem, a convex relaxed model has been proposed in [START_REF] Yıldızoğlu | A convex formulation for global histogram based binary segmentation[END_REF] and [START_REF] Swoboda | Variational image segmentation and cosegmentation with the wasserstein distance[END_REF] by considering the minimization over the convex set u ∈ [0, 1] N and the comparison of the non normalized color histograms instead of normalized distributions:

2)

The distribution µ 1 (resp. µ 0 ) is here normalized with respect to the mass of the histogram Hu (resp.

Once a global minimizer u * of problem (6.2) has been computed, a threshold ∈ [0; 1[ is required to get a binary segmentation of the image into two regions u (x) = if u * (x) ≥ and 0 otherwise. Notice that contrary to local data terms [START_REF] Chambolle | An algorithm for mean curvature motion. Interfaces and Free Boundaries[END_REF][START_REF] Nikolova | Algorithms for finding global minimizers of image segmentation and denoising models[END_REF], the thresholded solution may lead to a greater value of the functional, i.e. J(u (x)) ≥ J(u * ), since the coarea formula does not hold for general nonlinear functions D involving linear operators H, µ 1 1 T and µ 0 1 T . One other solution could be to rely on integer linear programming with methods such as branch-and cut [START_REF] Nemhauser | Integer and Combinatorial Optimization[END_REF], in order to directly estimate a global binary solution of problem 6.2.

In the following, the use of distances D related to Optimal Transport are described.

Optimal Transport between non normalized histograms

In order to have a robust distance between the non-normalized histogram, one would like to take D as the Wasserstein distance. Nevertheless, as the computation of Optimal Transport involves prohibitive computational costs, some approximations are now described.

Wasserstein-1 for 1 -D histograms

When considering 1 -D histograms µ and ν defined over a d = 1 dimensional space (for instance distributions of gray levels), the p Wasserstein distance reads [START_REF] Villani | Topics in Optimal Transportation[END_REF]:

Non-smooth Optimisation

We consider functionals that can be written as the sum of different terms being nonsmooth or non-convex, which is a common setup for problems in image processing. The objectif of this chapter is to review classical proximal splitting algorithms dedicated to the estimation of global or local minima of such functionals. Proximal splitting schemes are indeed first order optimization methods that present good performances for large scale problems, as the ones that will be considered in Parts II and III of this document. This chapter is inspired from [START_REF] Bauschke | Convex Analysis and Monotone Operator Theory in Hilbert Spaces[END_REF][START_REF]Introduction à l'optimisation convexe non lisse pour le traitement d'images[END_REF] in which rigorous proofs can be found.

A.1 Definitions

In the following, E is an Euclidean vectorial space of finite dimension equipped with an inner product •, • and a norm ||.|| = •, •

.

Definition 1 (Domain). Let f be a function defined from E to R = R ∪ +∞, we denote the domain of f as dom(f ) = {x ∈ E such that f (x) = +∞}.

Definition 2 (Convexity). Let f be a function defined from E to R, f is convex iff for all pairs (x, y)

Proposition 2. Let J be a convex, proper, l.s.c and coercive function defined on E with values in R, then J is bounded from below and admits at least one minimizer. If J is strictly convex then the minimizer is unique.

The FB algorithm can also be generalized to deal with a differentiable convex function and the sum of proximable convex ones [START_REF] Raguet | Generalized forward-backward splitting[END_REF].

When both functions f and g are non-smooth, the standard Douglas-Rachford (DR) algorithm [START_REF] Lions | Splitting algorithms for the sum of two nonlinear operators[END_REF][START_REF] Combettes | A Douglas-Rachford splitting approach to nonsmooth convex variational signal recovery[END_REF] relying on the proximal operators of the two functions can be considered.

Theorem 7 (Douglas-Rachford Algorithm). Let J = f + g be the sum of two convex, l.s.c, proper and coercive functions. Let µ ∈ [0, 2], τ > 0, x 0 ∈ E and (x n ) n∈N , (y n ) n∈N and (z n ) n∈N the sequences defined as

then there exists x * ∈ E minimizer of J such that the sequence (x n ) converges to x * .

The Alternating Direction Method of Multipliers (ADMM) is another classical optimization method for that purpose. This approach consists in introducing an auxiliary variable y ∈ E and dealing with the constraint y = x with a Lagrange multiplier z ∈ E. An augmented Lagrangian is then defined for τ > 0 and ∀(x, y, z)

One can show that problem (A.4) is equivalent to

Theorem 8 (ADMM Algorithm). Let J = f + g be the sum of two convex, l.s.c, proper and coercive functions. Let τ > 0, (y 0 , z 0 ) ∈ E × E and (x n ) n∈N , (y n ) n∈N et (z n ) n∈N the sequences defined as

then there exists z * ∈ E minimizer of J such that the sequence (z n ) converges to z * .

More details can be found in the reference papers [START_REF] Gabay | A dual algorithm for the solution of nonlinear variational problems via finite element approximation[END_REF][START_REF] Glowinski | Sur l'approximation, par éléments finis d'ordre un, et la résolution, par pénalisation-dualité d'une classe de problèmes de dirichlet non linéaires[END_REF]. The ADMM algorithm can be reformulated in term of proximal operators of f and g, which makes it equivalent to the DR algorithm [START_REF] Eckstein | On the Douglas-Rachford splitting method and the proximal point algorithm for maximal monotone operators[END_REF]. More precisely, the sequence of iterates of DR with µ = 1 for the primal probleml f + g are identical to the sequence of iterates of ADMM for the problem f * + g * [START_REF] Papadakis | Optimal transport with proximal splitting[END_REF].

A.2.2 Primal-Dual Proximal Algorithms

We now consider the minimization of two convex functions, one of them being composed with a linear operator K:

Let (x 0 , y 0 ) ∈ E × F and x0 = x 0 and the sequences

then the sequence (x n , y n ) converges to a solution (x * , y * ) of (A.14).

The time steps of this algorithm can be defined automatically with respect to the operator K by considering diagonal preconditioning matrices T and Σ instead of scalar values τ and σ [START_REF] Pock | Diagonal preconditioning for first order primal-dual algorithms in convex optimization[END_REF][START_REF] Lorenz | An inertial forward-backward algorithm for monotone inclusions[END_REF]. The algorithm can also deal with several functions including different linear operators i f i (K i x). Notice that better "worst case convergence rates" can be obtained by using adequate strategies involving extragradient steps [START_REF] Drori | A simple algorithm for a class of nonsmooth convex-concave saddle-point problems[END_REF], overrelaxation parameters [START_REF] Chambolle | On the ergodic convergence rates of a first-order primal-dual algorithm[END_REF] or inertial terms [START_REF] Lorenz | An inertial forward-backward algorithm for monotone inclusions[END_REF].

A final important point to underline is that the definition of the proximal operator is not limited to the L 2 norm (A.2) and any strongly convex function ψ can be considered:

Hence, depending of the function f , a good choice of the metric ψ can lead to an explicit expression of the proximal operator [START_REF] Chambolle | On the ergodic convergence rates of a first-order primal-dual algorithm[END_REF].

A.3 Non-convex functionals

We now return to the problem (A.4):

When f and or g are non-convex and one of the two functions is non-smooth, we would like to estimate a local minima of J.

A.3.1 Back to the Forward-Backward

In the non-convex setting, local minima can not be characterized by the subdifferential that it is only defined for convex functions. The set of critical points of J, that include its local minima, can nevertheless be characterized by some local properties derived from the Fréchet subdifferential. When J satisfies the so-called Kurdika-Lojasiewicz (KL) property, the authors of [START_REF] Attouch | Convergence of descent methods for semialgebraic and tame problems: proximal algorithms, forward-backward splitting, and regularized gauss-seidel methods[END_REF] have shown that the Forward-Backward Algorithm (A.5) converges to a critical point of J. We refer to [START_REF] Attouch | Convergence of descent methods for semialgebraic and tame problems: proximal algorithms, forward-backward splitting, and regularized gauss-seidel methods[END_REF][START_REF]Introduction à l'optimisation convexe non lisse pour le traitement d'images[END_REF] for more details on this technical assumption. Since the KL property is checked by all semi algebraic function, it it will be verified by all the considered image processing functionals introduced in the present document.

An acceleration of the Forward-Backward Algorithm has been proposed for this nonconvex setting [START_REF] Ochs | ipiano: Inertial proximal algorithm for nonconvex optimization[END_REF]. This approach relies on an inertial force, derived from the Nesterov acceleration dor the gradient descent [START_REF] Nesterov | Interior Point Polynomial Methods in Convex Programming: Theory and Algorithms[END_REF]: