Damiano Lombardi
email: damiano.lombardi@inria.fr

Emmanuel Maitre
email: emmanuel.maitre@imag.fr

Eulerian models and algorithms for unbalanced optimal transport

Keywords:

 to overcome this caveat are reviewed, and discussed. Due to the mix between transport and L 2 interpolation, these models can produce instantaneous motion at finite range. In this paper we propose new methods, parameter-free, for interpolating unbalanced densities. One of our motivations is the application to interpolation of growing tumor images.

Introduction

Context in image processing

Optimal transportation finds many applications, ranging from optimal design, to the description of geostrophic flows. One possible application is the interpolation and the registration of medical images. Many advanced registration techniques have been developed, based on the setting up of diffeomorphic maps, like the Large Deformation Diffeomorphic Metric Mapping (LDDMM), see for instance [START_REF] Miller | Geodesic shooting for computational anatomy[END_REF], or methods based on diffeomorphic demons, see [START_REF] Vercauteren | Diffeomorphic demons: Efficient nonparametric image registration[END_REF]. The method based on optimal transportation provides an interpolation that minimizes the kinetic energy of the map. This is not always realistic, depending upon the application, but one clear advantage of optimal transportation is that the solution exists and it is unique, and since pioneering works of Benamou and Brenier [START_REF] Benamou | A computational fluid mechanics solution of the Monge-Katorovich mass transfer problem[END_REF] it can be formulated as the optimization of a convex functional.

By structure, optimal transportation requires initial and final densities to be balanced, that is, of equal mass.

In real applications, this is seldom the case, either because images are projections of a 3D reality, or picture a growing object (e.g. tumor), or simply due to noise. Benamou [START_REF] Benamou | Numerical resolution of an "unbalanced" mass transfer problem[END_REF] proposed a way to tackle this latter problem. His approach was to find a compromise between an L 2 projection and an optimal transport. He considered the optimal transportation between the initial density and a final density of same mass, which was computed by minimizing its L 2 distance with the real final density. This produces a transport between the larger (in the L 2 meaning) part of the densities. Note that a weight is present in the method (see below) in front of this L 2 distance, and it is acknowledged in [START_REF] Benamou | Numerical resolution of an "unbalanced" mass transfer problem[END_REF] that this parameter choice could be problematic in some cases.

At a late stage of redaction of this article, we got aware of the work of Piccoli and Rossi [START_REF] Piccoli | Generalized Wasserstein distance and its application to transport equations with source[END_REF][START_REF] Piccoli | A generalized Benamou-Brenier formula for mass-varying densities[END_REF], who studied theoretically optimal transportation with source terms, generalizing Benamou's approach. In particular, they obtained an interpolated distance between the L 2 and Wasserstein metrics, and exhibit interesting properties of this new distance. Another approach to unbalanced optimal transport was proposed by Figalli and Gigli in [START_REF] Figalli | A new transportation distance between non-negative measures, with applications to gradients flows with Dirichlet boundary conditions[END_REF] (see also [START_REF] Ambrosio | A user's guide to optimal transport[END_REF]). Their approach is to consider the boundary as a source/sink of mass.

Our main concern is to provide a notion of generalized optimal transport interpolation between measures/densities of different mass that could apply to the study of tumors growth. Therefore we do not want to involve the picture boundary as in [START_REF] Figalli | A new transportation distance between non-negative measures, with applications to gradients flows with Dirichlet boundary conditions[END_REF], as far as mass is concerned. Rather, we will consider in our test cases isolated tumors (for real pictures) or isolated gaussians (for synthetic ones), between which growth occurs, but not by mass coming from the boundary. Likewise, we would like to avoid the "infinite speed" of L 2 interpolation that is inherent to the approach of [START_REF] Benamou | Numerical resolution of an "unbalanced" mass transfer problem[END_REF] and [START_REF] Piccoli | Generalized Wasserstein distance and its application to transport equations with source[END_REF][START_REF] Piccoli | A generalized Benamou-Brenier formula for mass-varying densities[END_REF]. In these two approaches, indeed, due to the fact that the resulting generalized interpolation mixes Wasserstein and L 2 metrics, some mass at positive distance of the initial tumor support could instantaneously appear on the interpolating path, which is undesirable for a tumor growth model.

In the following, we introduce several models that transport unbalanced densities. Our aim is to address the optimal transportation problem between two densities (images) which represent an object which has grown between two instants. These models are based on a modification of the projection method hidden in Benamou-Brenier algorithm, where a source term is added.

In the following sections, we first recall some basic facts about optimal transportation, in the balanced case. Then, we recall the modification for unbalanced densities proposed by Benamou and Piccoli-Rossi. We prove (prop. 1) that Piccoli-Rossi's distance is not properly speaking a generalization of Wasserstein distance since it does not recover genuine optimal transport in case of balanced densities. Next, we present our approach and discuss the general form of a source term, and give some explicit solution in simple (but relevant for applications) cases to illustrate their behavior (propositions 2, 3, 4). Next we present a simple modification of Benamou-Brenier algorithm (in 4.2) to handle a source term which does not depend on the optimization variables ρ, m. We then consider several source terms models: the first is built on the idea of reproducing an exponential growth, while the second uses a diffusion process to drive the motion of one density to the other. For both models we prove existence and uniqueness of minimizers. At last, several numerical test are provided. First we show on synthetic test cases that our algorithm performs well. The we apply it to real images of lung tumor growth, for which we find good agreement between our interpolated growth rate using the exponential model and experimental one.

Quick introduction to optimal transportation

Let Ω be an open bounded domain and let us consider the Monge problem of pushing one measure µ to another measure ν, through a transportation map which minimizes some cost. The standard setting assume that the measures µ and ν are absolutely continuous with respect to the Lebesgue measure, of densities ρ 0 and ρ 1 , nonnegative on Ω, and of equal mass:

Ω ρ 0 (x)dx = Ω ρ 1 (x)dx = 1.
In application to image processing, these densities will correspond to gray levels, and in general this condition would not be satisfied. A map T : Ω → Ω is a transfer map from ρ 0 to ρ 1 if for every subset A ⊂ Ω,

A ρ 1 (x)dx = {T (x)∈A} ρ 0 (x)dx. ( 1
)
If T is a C 1 mapping, then by a change of variables this is equivalent to

det(∇T (x))ρ 1 (T (x)) = ρ 0 (x),
which is under-determined. Let Γ(ρ 0 , ρ 1 ) be the set of mappings T transfering ρ 0 on ρ 1 . The L p Kantorovich-Wassertein distance between ρ 0 and ρ 1 is then defined by

d p (ρ 0 , ρ 1 ) p = inf T ∈Γ(ρ0,ρ1) |T (x) -x| p ρ 0 (x)dx.
The L p Monge-Kantorovitch problem (MKP) corresponds to find a mapping T such that this infimum is achieved.

In the case p = 2, the problem admits an unique solution (see e.g. Villani [START_REF] Villani | Topics in optimal transportation[END_REF] page 66), which is the gradient of a convex fonctional from Ω to R:

T (x) = ∇Ψ(x).
The convex function Ψ is solution of Monge-Ampère equation: This equation being highly nonlinear, numerical methods to solve the MKP problem based on discretization of the Monge-Ampère equation have already been investigated [START_REF] Loeper | Numerical solution of the Monge-Ampre equation by a Newton's algorithm[END_REF][START_REF] Dean | Numerical methods for fully nonlinear elliptic equations of the Monge-Ampère type[END_REF][START_REF] Benamou | Numerical solution of the second boundary value problem for the Elliptic Monge-Ampère equation[END_REF][START_REF] Benamou | Numerical Solution of the Optimal Transportation Problem via viscosity solutions for the Monge-Ampère equation[END_REF]. In application to image morphing problem, it is relevant to seek a time-dependent family of mappings T (•, t) transfering continuously ρ 0 to ρ 1 . In [START_REF] Benamou | A computational fluid mechanics solution of the Monge-Katorovich mass transfer problem[END_REF] the authors introduced a fluid mechanics formulation of MKP, by adding a new dimension to the original problem (the time). The idea is to consider an arbitrary time interval [0, 1] and all functions ρ(x, t) ≥ 0 and vector fields v(x, t) ∈ R n solution of the continuity conditions with prescribed initial and final densities:

det(D 2 Ψ)ρ 1 (∇Ψ(x)) = ρ 0 (x
∂ t ρ + div(ρv) = 0, ρ(x, 0) = ρ 0 (x), ρ(x, 1) = ρ 1 (x), (2)
and homogeneous Dirichlet conditions on ∂Ω. Then we have :

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

d 2 (ρ 0 , ρ 1 ) 2 = inf Ω 1 0 ρ(x, t)|v(x, t)| 2 dxdt
the infimum being taken on ρ, v verifying (2) in suitable functional spaces, see [START_REF] Villani | Topics in optimal transportation[END_REF].

Numerically, a first step is to make the problem convex by the change of variables (ρ, m) = (ρ, ρv). Then to solve a saddle-point problem based on a augmented Lagrangian method. In the last few years, others applications of optimal transportation methods to image analysis have been proposed. For instance, instead of solving the saddle-point problem directly, Angenent et al. derived a novel gradient flow for the computation of the optimal transport map [START_REF] Agenent | Minimizing flows for the Monge-Kantorovich problem[END_REF]. Unfortunately, all these methods require that the initial and final densities have the same mass. This can be seen directly on the mass conservation constraint (2) upon integration in time and space, using periodic or Dirichlet boundary conditions on v on ∂Ω.

The remaining of this paper is organized as follows: in the next section, after having recalled some existing solution to this problem of unbalanced densities, we consider a new optimal transport dealing with the different mass of densities. Several source terms are considered and compared, and numerical tests are performed.

Some existing models of unbanlanced mass transportation

First extensions of Kantorovitch norms to unbalanced measures was well described in [START_REF] Guittet | Extended Kantorovich norms: a tool for optimization[END_REF]: Kantorovich and Rubinstein [START_REF] Kantorovich | On a space of completely additive functions[END_REF], followed by a slight different definition of Hanin [START_REF] Hanin | An extension of the Kantorovich norm[END_REF] proposed generalized distances. However these true generalizations were mixing displacement distances and total variation, while we were looking for a pure displacement behavior.

The starting point of this article is the Benamou-Brenier fluid mechanics formulation of the L 2 -MKP [START_REF] Benamou | A computational fluid mechanics solution of the Monge-Katorovich mass transfer problem[END_REF]. Consider Ω = (0, 1) 2 with Dirichlet boundary conditions and a time interval [0, 1], we set Q m = Ω × (0, 1). In order to minimize the energy under the constraint (2), we first introduce the new variables ρ and m = ρv (into which the constraint expresses linearly) and we consider the (now convex) problem:

inf (ρ,m)∈C(ρ0,ρ1) Ω 1 0 |m| 2 2ρ dxdt
where

C(ρ 0 , ρ 1 ) = {(ρ, v), ∂ t ρ + div m = 0, ρ(•, 0) = ρ 0 , ρ(•, 1) = ρ 1 on Ω, m • n = 0 on ∂Ω × (0, 1)} . (3)
Note that upon space integration of the conservation equation in (3), we still get

Ω ρ 0 (x)dx = Ω ρ 1 (x)dx. (4)
In order to deal with unbalanced densities, Benamou proposed to somehow mix the L 2 and Wasserstein distances. Given a parameter γ > 0, one minimizes

d W (ρ 0 , ρ1 ) 2 + γ ρ1 -ρ 1 2 L 2
among all densities ρ1 of same mass as ρ 0 . While tests performed by Benamou showed that this algorithm allows to correctly compute interpolation between two densities with underlying noise, we observe that it is easy to imsart-generic ver. 2012/08/31 file: dlemOrth.tex date: October 7, 2014 find examples where, for γ large enough, ρ1 could be negative. This is the case for instance if ρ 0 is a gaussian of weight 1, while ρ 1 is the sum of two gaussians functions of weights 1 and 2 (we insist that this kind of examples was not under the scope of the method developed in [START_REF] Benamou | Numerical resolution of an "unbalanced" mass transfer problem[END_REF]). Then taking large γ would lead to lower significantly ρ 1 so that its smaller part could become negative. This is highly undesirable in the context we are considering.

A very related model of optimal transport with source term has been recently introduced by Piccoli and Rossi in two papers [START_REF] Piccoli | Generalized Wasserstein distance and its application to transport equations with source[END_REF][START_REF] Piccoli | A generalized Benamou-Brenier formula for mass-varying densities[END_REF]. In the first work, they provide a link between a transport equation with source term and a generalized Wasserstein distance, whereas the second paper introduces and studies the Benamou-Brenier formula in the case of unbalanced mass densities. More precisely, they proved that the generalized Wasserstein distance defined by

W a,b 2 (ρ 0 , ρ 1 ) 2 = inf ρ0, ρ1∈L 1 + (Ω), ρ0 1 = ρ1 1 a 2 ( ρ0 -ρ 0 1 + ρ1 -ρ 1 1 ) 2 + b 2 W 2 (ρ 0 , ρ1 ) 2
coincides with the generalized Benamou-Brenier formula, ie one has also:

W a,b 2 (ρ 0 , ρ 1 ) 2 = inf (ρ,v,S)∈C(ρ0,ρ1) a 2 1 0 S(•, t) 2 1 dt + b 2 1 0 Ω ρv 2 dxdt
where C(ρ 0 , ρ 1 ) = {(ρ, v, S), ∂ t ρ + div(ρv) = S, ρ(•, 0) = ρ 0 , ρ(•, 1) = ρ 1 } . Assumptions on the data are essentially identical to the no source term, see eg [START_REF] Villani | Topics in optimal transportation[END_REF]. They also prove that support localization result for S:

t∈[0,1] supp S(•, t) ⊂ t∈[0,1] supp ρ(•, t).
We remark that this results does not imply finite speed of propagation, that is, as this generalized distance is an interpolation between L 1 and W 2 , it is easy to construct examples of initial and final densities for which the interpolation ρ(•, t) has for t > 0 arbitrarily small a support at a fixed distance from the support of ρ 0 . Moreover, for balanced initial and final condition, the proposed interpolation does not in general recover the genuine Wasserstein distance, as illustrated in the following example in space dimension 1.

Proposition 1. Let Ω = (0, 1), and α ∈ (0, 1 2 ) such that 4a 2 α 2 < b 2 (1α). Consider the following initial and final balanced densities: ρ 0 (x) = 1 on (0, α) for 0 < α < 1 2 , and 0 elsewhere, and

ρ 1 (x) = ρ 0 (1 -x). Then W a,b 2 (ρ 0 , ρ 1 ) < b 2 d 2 (ρ 0 , ρ 1 )
and the optimal path for W a,b 2 differs from the Wasserstein interpolation.

Proof. Indeed the pure Wasserstein distance is 1α. Consider now the pointwise interpolation given by v = 0, ρ(x, t) = (1t)ρ 0 (x) + tρ 1 (x). This corresponds to S(x, t) = ρ 1 (x)ρ 0 (x). Thus

a 2 1 0 h(•, t) 2 1 dt + b 2 1 0 Ω ρv 2 dxdt = 4a 2 α 2 .
Therefore for a, b such that 4a 2 α 2 < b 2 (1-α), the Piccoli-Rossi generalized distance does not give the Wasserstein interpolation for mass balanced densities.

This caveat could somehow bring unphysical results when applied to real images, and one aim of this paper is to provide a generalized distance which recovers the genuine one for balanced densities.

At last, let us mention the user's guide to optimal transportation by Ambrosio-Gigli [START_REF] Ambrosio | A user's guide to optimal transport[END_REF], where the authors present a mass-varying optimal transport initially considered by Figalli and Gigli [START_REF] Figalli | A new transportation distance between non-negative measures, with applications to gradients flows with Dirichlet boundary conditions[END_REF]. Their aim is the following: knowing that the genuine Wasserstein metrics allows to define a solution to the heat equation as the flow of some energy with respect to that metric, the constant mass assumption leads to a Neumann boundary condition for this PDE. A natural question is to wonder how to modify the distance so that the resulting flow is a solution to an heat equation with Dirichlet boundary conditions. This is performed by restricting the transport condition on the interior of the domain, while leaving its boundary without condition. Existence of such a transport plan is proved, as well as properties of the resulting distance. This approach is not a remedy for our application to tumor growth, as far as they are isolated spots on a scanner picture. Indeed in that case, it would be quite unlikely that the growth occurs from a source coming from the boundary.

Models of unbalanced mass transport

General considerations

In this section a generic formulation of unbalanced mass transport is presented and some basic properties are investigated. In what follows, the Eulerian formulation of the optimal mass transport is adopted, that reads:

inf (ρ,m)∈C 1 0 Ω |m| 2 2ρ dxdt , (5)
C = {(ρ, m) | ∂ t ρ + ∇ • m = 0, ρ(x, 0) = ρ 0 , ρ(x, 1) = ρ 1 } . (6)
One of the basic properties of the optimal transport solution is the invariance with respect to time reflection, i.e. ρ(x, 1t), m(x, 1t) are solution of the problem when ρ 1 is transported in ρ 0 . This may be shown by simply considering the following transformation:

x ′ = x, (7)
t ′ = 1 -t, (8)
ρ ′ = ρ, (9)
m ′ = -m, (10)
λ ′ = -λ. (11)
The objective functional to be minimized is rewritten by perfoming this change of coordinates:

L ′ = 1 0 Ω |m ′ | 2 2ρ ′ + λ ′ (∂ t ′ ρ ′ + ∇ x ′ • m ′ ) dx ′ dt, (12)
that is formally equivalent to the original one. The change of coordinate leaves the functional (and the associated Euler-Lagrange equations) unchanged, so that the solution of the optimal transport will be invariant with respect to this transformation. Let us consider a generic source term, i.e. the constraint will be no longer homogeneous; instead it can be written as:

∂ t ρ + div m = S(x, t; ρ, m; ρ 0 , ρ 1 ), (13)
where S accounts for the mass variation and it may be a function of x, t as well as the variables (ρ, m) and all their derivatives in space and time, the initial and final density (for instance, for normalization purposes). Among all the possible source terms, it is meaningful to look for those preserving the symmetry property in time that characterizes the classical optimal transport problem. This is done by asking:

S(x, t, ρ, m; ρ 0 , ρ 1 ) = -S(x, 1 -t, ρ, -m; ρ 1 , ρ 0 ). (14)
Indeed, it may be checked that this condition is sufficient to leave the functional unchanged, so that the same argument shown for the balanced case may be adopted. Observe that on the right hand side the initial density is ρ 1 and the final one ρ 0 . As well, a natural condition on S would be to vanish when ρ 0 and ρ 1 have the same mass, so that we recover classical optimal transportation :

Ω ρ 0 dx = Ω ρ 1 dx ⇒ {S(x, t, ρ, m; ρ 0 , ρ 1 ) = 0} . ( 15
)
The Euler-Lagrange equations for the objective functional to be minimized read:

m ρ -∇λ -λ δS δm = 0, (16)
∂ t λ + |m| 2 2ρ 2 + λ δS δρ = 0, (17)
imsart-generic ver. 2012/08/31 file: dlemOrth.tex date: October 7, 2014

The algorithm presented in this work approximates the solution of Eq.( 16)-( 17) for different source models, of increasing complexity. The affine case amounts to consider a non homogeneous mass conservation, with a constant prescribed source term. While this case is of limited interest for application to real images, we will show that the Benamou-Brenier algorithm easily adapts to that situation. Moreover, for less trivial and time dependent source terms, we will use this algorithm, with an explicit scheme (i.e. by taking the source term at the previous time step).

The first non constant source term we will consider is an exponential model, which is is interesting for modeling the corresponding behavior of tumor growth. We will show explicit solutions for simple cases where the initial and final densities are linked either by a translation/scaling or a affine transformation/scaling.

The second source term considered will somehow use a dual Sobolev norm to estimate the distance between the density pair, by solving a stationary Laplace equation, and use the corresponding flux to drive the mass growth. For that model we will be able to show existence of solution by rephrasing the model as a Monge problem on a manifold and using results of McCann [START_REF] Mccann | Polar factorization of maps on Riemannian manifolds[END_REF].

At last, we will consider a normal growth model where the source term is proportional to the modulus of ∇ρ. While we will not be able to theoretically prove existence on that model, it will turn to be the best choice for the modeling of tumor growth in some situations.

Remarks on the qualitative behavior of the solutions

In this section some qualitative remarks on solutions are proposed. First, as constant speed transport is a particular solution to genuine optimal transport [START_REF] Benamou | A computational fluid mechanics solution of the Monge-Katorovich mass transfer problem[END_REF], we found it interesting to investigate whether the solution of an unbalanced optimal transport may be merely a translation at constant speed and a rescaling. This is answered by the following proposition.

Proposition 2. Let ρ(x, t) = ρ 0 (x -at)µ(t)
, where a is a constant vector field (so that ∂ t a = 0 and ∇a = 0) and µ(t) a smooth scaling factor. If (ρ, m) with m = ρa verifies ( 13) and the Euler-Lagrange equations [START_REF] Mccann | Polar factorization of maps on Riemannian manifolds[END_REF][START_REF] Peyre | Non-asymptotic equivalence between W 2 distance and H -1 norm[END_REF], then a = 0 or µ is constant in time.

Proof. In order for the mass conservation equation to be satisfied it is necessary that:

S = ρ∂ t log µ(t). (18)
From ( 16) we get m ρ = ∇λ which plugged in [START_REF] Peyre | Non-asymptotic equivalence between W 2 distance and H -1 norm[END_REF] gives:

∂ t λ + |∇λ| 2 2 + λ∂ t log µ = 0. (19)
As m = ρa we have a = ∇λ; hence, by taking the gradient of Eq.( 19), which is an Hamilton-Jacobi equation, the following condition is found:

a∂ t log µ = 0, (20)
that implies that a = 0 or µ = const.

This means that a pure translation is a solution of a genuine optimal transport problem, a pure scaling is a solution of a differential equation, but a combination of the two can not be solution of an optimal transportation with some source term S. The next paragraph will provide an analytical solution where the translation vector is time varying.

Translation and scaling for an exponential model of growth

Let us look for a particular analytic solution in the case of an exponential source unbalanced optimal transport. Proposition 3. Consider the problem of finding an optimal plan between the densities:

ρ(x, 0) = ρ 0 , ρ(x, 1) = ρ 0 (x -a)e c , (21)
imsart-generic ver. 2012/08/31 file: dlemOrth.tex date: October 7, 2014

where ρ 0 is a positive smooth function, a is a constant vector field and c ∈ R. When an exponential model source given by S = cρ is considered, then

ρ(x, t) = ρ 0 x - 1 -e -ct 1 -e -c a e ct , m = ρa, (22)
is verifying the mass conservation ( 13) and the Euler-Lagrange equations [START_REF] Mccann | Polar factorization of maps on Riemannian manifolds[END_REF][START_REF] Peyre | Non-asymptotic equivalence between W 2 distance and H -1 norm[END_REF].

Proof. For this source term, ( 13) and [START_REF] Mccann | Polar factorization of maps on Riemannian manifolds[END_REF][START_REF] Peyre | Non-asymptotic equivalence between W 2 distance and H -1 norm[END_REF] are:

∂ t ρ + ∇ • m = cρ, (23) m ρ
-∇λ = 0 (24)

∂ t λ + |∇λ| 2 2 = -cλ. ( 25
)
By setting m = ρv, we get v = ∇λ. By taking the gradient of the equation for the lagrangian multiplier, this velocity should verifiy:

∂ t v + v∇v = -cv. ( 26
)
Introducing the flow of v, that is (ξ, t) → X(ξ, t) solutions of:

∂ t X = v(X(ξ, t), t), X(ξ, 0) = ξ, (27)
we get by plugging v(X(ξ, t), t) in (26) that

v(X(ξ, t), t) = v 0 (ξ)e -ct
for some initial velocity v 0 . Therefore there exist a constant in time d(ξ) such that

X(ξ, t) = - 1 c v 0 (ξ)e -ct + d(ξ).
We already have X(ξ, 0) = ξ, and a natural condition is to look for a solution such that X(ξ, 1) = ξ + a. Then we find that

v 0 (ξ) = ac 1 -e -c and X(ξ, t) = ξ + 1 -e -ct 1 -e -c a. ( 28
)
Let us observe that v obtained from X by ( 27) is constant in space. The continuity equation may thus be re-written as:

∂ t ρ + v • ∇ρ = cρ, (29)
and by introducing σ = ρe ct , this equation reduces to an homogeneous transport equation. The general solution is therefore written as:

ρ(x, t) = ρ 0 x - 1 -e -ct 1 -e -c a e ct . ( 30
)
Let us remark that in the limit of a vanishing c, that corresponds to a classical optimal transport, a pure translation is recovered:

lim c→0 ρ(x, t) = ρ 0 (x -at). (31)

Affine mapping: unbalanced homothety

In that section we consider an initial density function ρ 0 and a final density function obtained from ρ 1 by an affine mapping and a dilatation.

ρ(x, 0) = ρ 0 (x), ρ(x, 1) = µρ 0 (M x), M ∈ R d×d , µ ∈ R + , ( 32
)
where M is an invertible matrix whose entries do not depend upon space and time coordinates, and d is the space dimension. When an exponential source term is considered, namely S = cρ with c = log(µ/ det(M )), then

ρ(x, t) = e ct ρ 0 (Bx) det(B) (33)
is solution to the mass conservation ( 13) and the Euler-Lagrange equations [START_REF] Mccann | Polar factorization of maps on Riemannian manifolds[END_REF][START_REF] Peyre | Non-asymptotic equivalence between W 2 distance and H -1 norm[END_REF], with

B(t) := I + 1 -e -ct 1 -e -c A -1
and

A = M -1 -I. ( 34
)
Proof. For the particular case of an exponential model we have:

∂ t ρ + ∇ • (ρv) = cρ, (35)
∂ t λ + |∇λ| 2 2 = -cλ. ( 36
)
After integration on the space domain of the mass conservation equation, the following relation is obtained:

D t Ω ρ dΩ = c Ω ρ dΩ, (37)
Integrating in time between t = 0 and t = 1, this provides:

Ω µρ 0 (M x) dx = e c Ω ρ 0 dx. (38)
By performing a change of coordinate in the integral on the left hand side, and by using the fact that the matrix M is given and constant in space, the exponent c may be expressed as function of the determinant of the matrix and the scaling factor µ:

c = ln µ det(M ) . ( 39
)
As in the former proposition, the characteristic lines

of v = m ρ verify X(ξ, t) = -1 c v 0 (ξ)e -ct + d(ξ) which with X(ξ, 0) = ξ gives X(ξ, t) = ξ + 1 c v 0 (ξ)(1 -e -ct ).
To reach the final density we look for a solution such that X(ξ, 1) = M -1 ξ, which gives

v 0 (ξ) = c M -1 ξ -ξ 1 -e -c , X(ξ, t) = ξ + (M -1 ξ -ξ) 1 -e -ct 1 -e -c = I + 1 -e -ct 1 -e -c A ξ (40)
where we set (I + A) -1 = M . Therefore the backward characteristics are given by

Y (x, t) = B(t)x where B(t) := I + 1 -e -ct 1 -e -c A -1
.

Standard result on conservation laws state that

∂ t [ρ(X(ξ, t), t) det ∇X(ξ, t)] = (∂ t ρ + div(ρv))(X(ξ, t), t)
which turn equation (35) to a first order ODE which is trivially integrated. Using Y (ξ, t) = X -1 (ξ, t) we obtain

ρ(x, t) = e ct ρ 0 (Y ) det(∇ x Y ), (41)
that, after substitution of the expression of Y becomes:

ρ(x, t) = e ct ρ 0 (B(t)x) det(B(t)), (42)
Remark that as B(0) = I and B(1) = M , the constraints on the initial and final densities are recovered.

In the following sections some particular examples of source terms are detailed, that lead to different solutions for the unbalanced optimal transport. Their properties will be investigated by means of numerical experiments.

imsart-generic ver. 2012/08/31 file: dlemOrth.tex date: October 7, 2014 4. Variables independent source: affine constraint

Formulation

We consider the case where we do not impose (4) anymore. This is particularily relevant in the case where we are interpolating between two images of a growing tumor, for instance. A natural idea is to add a source term in the mass conservation constraint. As this mass growth is supposed to hold on the boundary of domains of homogeneous densities which represents structures in the image, a first guess would be to consider a mass conservation constraint where the source term S does not depend on ρ and m. Then the constraint on mass conservation is modified like:

∂ t ρ + div(m) = α (43)
where α : Ω × (0, 1) → R is a given function verifying

1 0 Ω αdx = Ω ρ 1 (x) -ρ 0 (x)dx.
The simplest form could be to consider is a constant α, which gives:

α = 1 |Ω| Ω ρ 1 (x) -ρ 0 (x)dx (44)
but while it verifies ( 14)-( 15) we will see that this is not the best choice for applications (due to the fact that it is not localized), and therefore we keep a possibly time and space dependant α in the following. The associated Lagrangian is given by

L(φ, ρ, m) = 1 0 Ω |m| 2 2ρ -ρ∂ t φ -m • ∇φ -αφdxdt - Ω φ(0, x)ρ 0 (x) -φ(1, x)ρ 1 (x)dx. (45)
Given two densities ρ 0 et ρ 1 , the minimization problem is equivalent to the saddle-point problem:

inf (ρ,m)∈Cα(ρ0,ρ1) sup φ L(φ, ρ, m),
where

C α (ρ 0 , ρ 1 ) = {(ρ, v), ∂ t ρ + div m = α, ρ(•, 0) = ρ 0 , ρ(•, 1) = ρ 1 on Ω, m • n = 0 on ∂Ω × (0, 1)} . ( 46
)
Arguing as in [START_REF] Benamou | A computational fluid mechanics solution of the Monge-Katorovich mass transfer problem[END_REF] we introduce dual variables (a, b) ∈ R × R d such that

|m| 2 2ρ = sup (a,b)∈K a(t, x)ρ(t, x) + b(t, x) • m(t, x), with K = (a, b) : R × R 2 → R × R 2 , a + 1 2 |b| 2 ≤ 0 on R × R 2 .
For sake of clarity, we set µ = (ρ, m) and q = (a, b), and introduce the support function of K, F such that F (q) = 0 for q ∈ K and F (q) = +∞ otherwise. Therefore we have

|m| 2 2ρ = sup q∈K µ • q = sup q -F (q) + µ • q. At last we set G(φ) = Ω φ(0, x)ρ 0 (x) -φ(1, x)ρ 1 (x)dx.
Still following [START_REF] Benamou | A computational fluid mechanics solution of the Monge-Katorovich mass transfer problem[END_REF], we show that our saddle point problem can be written as

sup µ inf φ,q F (q) + G(φ) + µ, ∇ t,x φ -q + α, φ , (47)
where the brackets stand for the L 2 (Q m ) scalar product, the variables µ, q are taken in L 2 (Q m ) d+1 , and φ in H 1 (Q m ). We now aim at finding a saddle-point of this problem which corresponds to a standard form of [START_REF] Fortin | Augmented Lagrangian Methods: Applications to the Solution of Boundary Value Problems[END_REF] in order to apply augmented Lagrangian techniques. The formal optimal condition for this problem are:

imsart-generic ver. 2012/08/31 file: dlemOrth.tex date: October 7, 2014

             ∂ t φ + |m| 2 2ρ 2 = 0 in [0, 1] × Ω ∂ t ρ + div m = α in [0, 1] × Ω m ρ = ∇φ in [0, 1] × Ω ρ(0, .) = ρ 0 in Ω ρ(1, .) = ρ 1 in Ω
Observing that the variable m can be eliminated, the optimality conditions can be rewritten in term of ρ, φ and c as:

       ∂ t φ + |∇φ| 2 2 = 0 in [0, 1] × Ω ∂ t ρ + div(m) = α in [0, 1] × Ω ρ(0, .) = ρ 0 in Ω ρ(1, .) = ρ 1 in Ω (48)
Therefore the optimal mass transfer still follows straight lines. We then define the augmented Lagrangian by introducing r > 0:

L r (φ, q, µ) = F (q) + G(φ) + µ, ∇ t,x φ -q + α, φ + r 2 ∇ t,x φ -q, ∇ t,x φ -q . ( 49
)

Algorithm

We consider the following iterative algorithm to compute this saddle point numerically: This algorithm builds from (φ n-1 , q n-1 , µ n , c n-1 ) the next iterate, and is very close to the original Benamou-Brenier algorithm. We just describe the differences in the three steps.

Step A φ n = arg min L r ( . , q n-1 , µ n ) This still amounts to solve a Poisson equation, but now with an extra contribution coming from α. Namely, taking the differential with respect to φ gives G(φ) + µ n , ∇ x,t φ + α, φ + r ∇ x,t φ nq n-1 , ∇ x,t φ = 0, ∀φ which, for Dirichlet boundary conditions gives, following [START_REF] Benamou | Numerical resolution of an "unbalanced" mass transfer problem[END_REF]:

-r∆ x,t φ n = div x,t (µ nrq n-1 )α (x, t) ∈ Ω×]0, 1[ with non homogeneous Neumann boundary conditions in space and time:

r∂ t φ n (x, 0) = ρ 0 (x) -ρ n (0, x) + ra n-1 (x, 0) r∂ t φ n (x, 1) = ρ 1 (x) -ρ n (1, x) + ra n-1 (1, x) (50)
r∂ n φ n (x, t) = rb n-1 (x) • n -m n • n on ∂Ω × (0, 1) (51)
Step B q n = arg min L r (φ n , . , µ n ) is identical to [START_REF] Benamou | Numerical resolution of an "unbalanced" mass transfer problem[END_REF], i.e. a pointwise projection on a paraboloid.

Step C µ n+1 = arg max L r (φ n , q n , . ) is identical to [START_REF] Benamou | Numerical resolution of an "unbalanced" mass transfer problem[END_REF].

Source proportional to a scalar field

In this section another source model is investigated of the form:

S(x, t; ρ, m; ρ 0 , ρ 1 ) = -ρ∂ t Γ -m • ∇Γ = -ρD t Γ, (52)
where Γ(x, t) is a given scalar field and D t is the total (lagrangian) derivative. Two different cases will be investigated: an exponential type of growth and a heat flux guided growth. The first will correspond to Γ(x, t) = ct for c ∈ R which will be determined from initial and final densities (section 5.1). The second will correspond to a time-independent Γ which will be built on an heat flux model (section 5.2). This seemingly strange expression of source term is however natural: first it is a quite general linear form in (ρ, m) (thus preserving convexity) and in the two mentioned cases (i.e. Γ constant in time or space) we are able to prove existence and uniqueness of solution to the corresponding optimization problem.

Remark 1. Mass sources could also be use to model obstacles. This could be useful in image interpolation with constraints : for instance when some part of the image is fixed. A mass source of expression S = -ρD t Γ, with Γ = Γ(x, t) may act as an obstacle. If the obstacle is fixed in space, then Γ is time independent, thus S = -m • ∇Γ. This may be deduced by performing a change of variable in the objective functional. In particular, let us consider the lagrangian:

L = 1 0 Ω |m| 2 2ρ + φ ∂ t ρ + ∇ • m + ρ DΓ Dt dx dt. (53)
A density variable σ defined as:

ρ := σe -Γ , (54)
is introduced, that inserted into the lagrangian transform the problem into:

L = 1 0 Ω 1 2 σe -Γ |v| 2 + λ (∂ t σ + ∇ • (σv)) dx dt, ( 55
)
where λ is the Lagrange multiplier associated to the constraint on σ. Let us remark that in this case the constraint is the usual one and the kinetic energy is modified by an isotropic metric term involving the exponential of the scalar field. The Euler-Lagrange equation for this system reads:

∂ t λ + v • ∇λ = 1 2 e -Γ v 2 , ( 56
) v = e Γ ∇λ, (57)
that, after substitution reduces to:

∂ t λ + 1 2 e Γ |∇λ| 2 = 0. (58)
In conclusion, every source term of the form ρD t Γ acts in a dual manner. It may be considered as a source term or a metric factor.

Exponential model of growth

Introduction

For this first model, let us consider a source of the form:

S = -cρ (59)
where c ∈ R is a constant. This corresponds to Γ(x, t) = ct in (52). Integrating ∂ t ρ + div m = -cρ in space gives, using the homogeneous boundary conditions on m,

d dt Ω ρdx = -c Ω ρdx. (60)
The expression for c may be computed a priori, depending only on the initial and the final mass only. Indeed upon integration of the first order ODE (60) we get:

c = log Ω ρ 0 dx Ω ρ 1 dx . ( 61
)
Note that the source term (59) with c given by (61) verifies ( 14)-( 15).

Remark 2. Another way to derive an expression for c is to integrate both members of (60), and this leads to a growth rate which is still constant in space and time, but depends nonlinearly on ρ:

c[ρ] = Ω ρ 0 -ρ 1 dx 1 0 Ω ρ dx . ( 62
)
Therefore, the first c can be computed a priori while the second would have to be reevaluated at each iteration. One could wonder whether the linear growth model and the nonlinear one would give the same optimal path, if it exists. We have equivalence of these two models of growth:

imsart-generic ver. 2012/08/31 file: dlemOrth.tex date: October 7, 2014 Proposition 5. Let ρ 0 , ρ 1 ∈ L 1 (Ω) be nonnegative, with positive integrals on Ω, and define the sets

C 1 = (ρ, m) ∈ L 1 ((0, 1) × Ω) × L 1 (0, 1; W 1,1 (Ω)), ∂ t ρ + div m = -cρ, c = log Ω ρ 0 dx Ω ρ 1 dx ρ(•, 0) = ρ 0 , ρ(•, 1) = ρ 1 , m • n = 0 on ∂Ω × (0, 1)} (63)
and

C 2 = (ρ, m) ∈ L 1 ((0, 1) × Ω) × L 1 (0, 1; W 1,1 (Ω)), ∂ t ρ + div m = -c[ρ]ρ, c[ρ] = Ω ρ 0 -ρ 1 dx 1 0 Ω ρ dx ρ(•, 0) = ρ 0 , ρ(•, 1) = ρ 1 , m • n = 0 on ∂Ω × (0, 1)} (64) Assume C 1 = ∅. Then C 1 = C 2 .
Proof. Let (ρ, m) ∈ C 1 . From the positiveness of the integrals of initial and final densities, integrating (60) from 0 to t and from t to 1 for t ∈ (0, 1) gives (assuming c ≤ 0 without loss of generality):

0 < Ω ρ 0 dx ≤ Ω ρ(x, t)dx ≤ Ω ρ 1 dx
on (0, 1).

Therefore we can compute

Existence and uniqueness of the solution

In this section the existence and uniqueness of the solution for the exponential source are investigated. The proof of the proposition is based on the following: first, the objective functional defining the transport problem is transformed. It is shown that the transport is equivalent, in the new variables, to an homogeneous transport with a time dependent metric. Then, the geodesics of the transport are studied, allowing to show that the distance between the densities is proportional to the classical Wasserstein distance. In particular, as done for the classical Benamou-Brenier transport, it is shown than the objective functional in eulerian form, which is larger or equal to the transport cost (by Jensen inequality) coincides with the transport cost for the geodesic velocity field. The result relies on the following proposition: Proposition 6. Let Ω ⊆ R d be a convex. The geodesic curve γ(t) : [0, 1] → Ω with respect to the metric tensor g ij := e -ct δ ij , c ∈ R is defined as:

γ(t) = arg min γ 1 0 e -ct | γ| 2 dt. (65)
We prove that γ(t) = a(c, t)γ(1) + b(c, t)γ(0) for a, b two real functions of (c, t), and that the geodesic distance squared

d 2 (γ(0), γ(1)) ∝ |γ(1) -γ(0)| 2 .
Proof. The proof is based on a direct computation. The equation minimizing the cost in Eq.( 65) reads:

d dt e -ct d dt γ = 0, (66)
whose integration between γ(0) and γ(1) provides:

γ(t) = e c e c -1 γ(0) - 1 e c -1 γ(1) + e ct e c -1 (γ(1) -γ(0)), (67)
which prove the first part of the proposition. Remark that the geodesics are straight line parametrized at non-constant speed. The time derivative of γ is computed and inserted into the expression of the distance, providing:

d(γ(0), γ(1)) 2 = c e c -1 |γ(1) -γ(0)| 2 , (68)
that concludes the proof.

imsart-generic ver. 2012/08/31 file: dlemOrth.tex date: October 7, 2014

The result of this proposition will be used to proof the existence and uniqueness result in the following:

Proposition 7. Let ρ 0 , ρ 1 be two given density distributions. The solution of the optimal transport problem with exponential mass source defined by Eq.( 59) exists and it is unique.

Proof. The objective functional defining the problem is:

L = 1 0 Ω 1 2 ρv 2 + λ(∂ t ρ + ∇ • (ρv) + cρ) dx dt. (69)
The following change of variable is performed, σ := ρe ct , leading to the equivalent problem:

L = 1 0 Ω 1 2 σe -ct v 2 + φ(∂ t σ + ∇ • (σv)) dx dt, (70)
which is an homogeneous transport, with a kinetic energy depending upon a time variable isotropic metric. For the Jensen inequality, it holds:

1 0 Ω0 1 2 σ 0 (ξ)e -ct v 2 dξ dt ≥ Ω0 1 2 σ 0 (ξ)d 2 (ξ, X(ξ)) dξ, (71)
where d 2 (ξ, X) is the geodesic distance squared between ξ ∈ Ω 0 and x = X(ξ) ∈ Ω 1 , the initial and final configurations. Let us prove that the equality holds for a velocity field v = Ẋ which is related to the solution of the optimal transportation problem between σ 0 and σ 1 . To do so, the result of Proposition 6 is used. The minimizer of the geodesic distance under the mass constraint satisfies:

X = arg min X∈C Ω0 1 2 σ 0 (ξ)d 2 (ξ, X(ξ)) dξ = arg min X∈C c e c -1 Ω0 1 2 σ 0 (ξ)|X -ξ| 2 dξ, (72)
C = {X|σ 0 (ξ) = σ 1 (X(ξ) det(∇ ξ X))} . ( 73
)
The solution is a classical optimal transport between σ 0 and σ 1 . There exists a unique minimizer and X(ξ, 1) = ξ + ∇ ξ Φ, where ξ 2 2 + Φ(ξ) is a convex function. This solution is parametrized at non-constant speed for the time mapping X(ξ, t) to be geodesic:

X(ξ, t) = ξ + e -ct -1 e c -1 ∇ ξ Φ ⇒ v = Ẋ = ce ct e c -1 ∇ ξ Φ. ( 74
)
The expression of the velocity is introduced into the Eq.( 71) and the following holds:

min Ẋ 1 0 Ω0 1 2 e -ct σ 0 Ẋ2 dξ dt = min X Ω0 1 2 σ 0 d 2 (ξ, X) dξ = 1 2 c e c -1 W 2 (σ 0 , σ 1 ), ( 75
)
where Ẋ is chosen as in Eq.( 74). Thus, the problem of minimizing the objective functional in Eq.( 69) admits a unique solution.

Remark 3. If ρ 0 = ρ 1 e c , the change of variable provides σ 0 = σ 1 ⇒ v = 0. This corresponds to a trivial solution, in which there is no transport, but ρ 0 , ρ 1 are obtained by a pointwise interpolation.

Heat flux guided model

The counterpart of the exponential model of growth consists in taking S = -m • ∇Γ(x), where Γ may be either a given field depending upon some known information associated to the problem, or a quantity to be determined as function of the problem data (i.e. ρ 0 , ρ 1 ) in order to set up a parametric free model. A perspective on the modeling of constraints has been presented at the beginning of this section. A preliminary constraint on Γ is derived.

Proposition 8. Let Γ ∈ H 1 (Ω). It can be an admissible source potential if ∆ρ, e Γ = 0, where , denotes the standard L 2 scalar product and ∆ρ = ρ 1ρ 0 .

Proof. Let us introduce σ(x, t) := ρ(x, t)e Γ(x) . Hence:

∂ t ρ + ∇ • ρv + m • ∇Γ = 0 ⇒ ∂ t σ + ∇ • (σv) = 0, (76)
so that σ satisfies an homogeneous constraint. This implies:

Ω σ(x, 0) dΩ = Ω σ(x, 1) dx ⇒ Ω (ρ 1 -ρ 0 )e Γ dx = 0, ( 77
)
the exponential of the source Γ have to be orthogonal to the density difference with respect to the L 2 scalar product Among all the possible admissible sources Γ, a heat-like solution is adopted. The source Γ is the solution of:

(Γ * , µ * ) = arg inf Γ sup µ Ω 1 2 |∇Γ| 2 - |δM | δM ∆ρΓ dx -µ Ω ∆ρe Γ dx, ( 78
)
where the scalar µ is the lagrange multiplier enforcing the orthogonality constraint. Homogeneous Dirichlet boundary conditions are used for Γ.

The equations minimizing the functional read:

-∇ 2 Γ = |δM | δM + µe Γ ∆ρ, in Ω, (79)
Γ = 0, on ∂Ω, (80)
Ω ∆ρe Γ dΩ = 0. (81)
Proposition 9. Let Γ be the solution of Eq.(81). Then, Γ satisfies the constraint expressed in the Proposition 8 and the solution of the optimal transport respects the time-reversal symmetry.

Proof. The constraint expressed in Proposition 8 is directly imposed in the third equation of the system Eq.( 81). When δM = 0, that is when the initial and final densities have the same mass, Γ = 0, identically, and a classical optimal transport is recovered. To verify that the time-reversal symmetry of optimal transport is satisfied, the following transformation is applied to the solution of the system:

∆ρ ′ = -∆ρ (82) µ ′ = -µ ′ . ( 83
)
From the first transformation, the sign of δM changes and the following is obtained:

-∇ 2 Γ ′ = - |δM | δM -µe Γ ′ (-∆ρ), in Ω, (84)
Γ ′ = 0, on ∂Ω, (85)
so that Γ ′ = Γ and the time reversal symmetry is preserved.

The boundary integral of the normal derivative of Γ is related to the absolute value of the mass difference between the densities. Indeed, the integral over the whole domain of the solution leads to:

- ∂Ω ∂ n Γ dS = |δM |. ( 86
)
Note that the source Γ encodes the distance between the supports of the densities (it is related to the H -1 distance between them) and the mass difference.

The problem of finding Γ solution of Eq.( 81) is solved by an Uzawa augmented lagrangian method, the starting value adopted for the lagrangian multiplier µ being µ = 0.

Existence and uniqueness of the solution

The existence and uniqueness of the solution of the optimal transport with a heat-flux mass source are investigated. The strategy of the proof is similar to that used for the exponential source. First, a change of variable is performed in such a way that the problem is transformed into an homogeneous transport with a metric. Then, a study of the regularity of the metric allows to apply directly a result of existence and uniqueness (see [START_REF] Mccann | Polar factorization of maps on Riemannian manifolds[END_REF] for the theorem).

Proposition 10. Let Γ be a strong solution of Eq.81. Then, the solution of:

(ρ, v, λ) = inf ρ,v sup λ 1 0 Ω 1 2 ρv 2 + λ(∂ t ρ + ∇ • (ρv) + ρv • ∇Γ) dx dt, (87)
exists and it is unique.

Proof. As seen in Eq.( 76) it is possible to rewrite this model through a change of variables in such a way that the constraint, in the new variables, is the classical homogeneous one. The objective functional, in the new variables, may be written as:

L = 1 0 Ω 1 2 σe -Γ v 2 + φ(∂ t σ + ∇ • (σv)) dx dt, (88)
which is an homogeneous transport with a scalar metric factor e -Γ . Let us remark that in the case of balanced densities Γ = 0 and the classical optimal transportation problem in Benamou-Brenier formulation is recovered.

The problem can be formulated as an optimal mass transport on a manifold. Indeed, let x = X(ξ, t)|X(ξ, 0) = ξ and v(X(ξ, t), t) = ∂ t X(ξ, t). Then the objective functional can be rewritten as:

L = 1 0 Ω0 1 2 σ(X)e -Γ(X) Ẋ2 + ψ(ξ)(σ 0 (ξ) -σ(X) det(∇ ξ X)) det(∇ ξ X) dξ dt. (89)
By applying the mass conservation equation in Lagrangian form, the time integral of the kinetic energy transforms into:

K = 1 0 Ω0 1 2 σ 0 (ξ)e -Γ(X) Ẋ2 dξ dt. ( 90
)
The problem is thus equivalent to the minimization of K under the homogeneous mass constraint between σ 0 and σ 1 . This case may be considered as an optimal transportation on a Riemannian manifold, whose metric tensor is g ij = e -Γ(X) δ ij . Let us restrict to the strong solutions of Eq.( 81), so that Γ ∈ C 2 (Ω). This is a sufficient condition for the manifold to be C 3 smooth, i.e. for the metric tensor to be twice continuously differentiable. Thus, by applying the results of [START_REF] Mccann | Polar factorization of maps on Riemannian manifolds[END_REF], the solution of this problem exists and it is unique. Moreover, in the minimum, the kinetic energy integral satisfies:

1 0 Ω0 1 2 σ 0 (ξ)e -Γ(X) Ẋ2 dξ dt = Ω0 1 2 σ 0 (ξ)d 2 (ξ, X) dξ, (91)
where d is the geodesic distance on the manifold and X is such that σ 0 (ξ) is mapped into σ 1 (X(ξ)) by the exponential map on the manifold (see [START_REF] Mccann | Polar factorization of maps on Riemannian manifolds[END_REF] for more details).

Source proportional to the modulus of the density gradient: nonlinear constraint

Another formulation for the source term would be to introduce some nonlinear growth term in α, which leads to a normal growth. We propose the following form:

α = β|∇ρ|, where β = Ω ρ 1 -ρ 0 dx 1 0 Ω |∇ρ|dx . ( 92
)
This choice is justified as follows: it merely says that mass variation is more localized on regions where ρ is varying. This seems natural to see a growing set as gaining mass on the boundary. The β term ensures that the gain of mass is compatible with the difference of mass between initial and final densities. The corresponding set of constraints

C α (ρ 0 , ρ 1 ) = (ρ, v), ∂ t ρ + div m = Ω ρ 1 -ρ 0 dx 1 0 Ω |∇ρ|dx |∇ρ|, ρ(•, 0) = ρ 0 , ρ(•, 1) = ρ 1 on Ω, m • n = 0 on ∂Ω × (0, 1)} . ( 93
)
seems not easy to deal with numerically. However, in the context of our augmented Lagrangian approach, we can use a quasi-static formulation which, at each iteration, amounts to use for α its expression in terms of ρ n . This means that our algorithm reduces to the original Benamou-Brenier algorithm where we just change step 1 to the Nonlinear step 1 φ n = arg min L r ( . , q n-1 , µ n ) which amounts to solve the Poisson equation,

-r∆ x,t φ n = div x,t (µ n -rq n-1 ) -Ω ρ 1 -ρ 0 dx 1 0 Ω |∇ρ n |dx |∇ρ n | (x, t) ∈ Ω×]0, 1[
with the same non homogeneous Neumann boundary conditions in space and time.

Numerical tests

In this section, we will provide some numerical study of the exponential and heat-flux guided source models.

On synthetic examples (gaussians of different masses in several configurations) and more realistic ones, we show that they behave quite well. We do not compare them with Benamou or Piccoli-Rossi's approaches. Indeed, Benamou's extension, as mentioned in section 2 was not developed to deal with big difference of mass between densities, but rather with noisy data. The other approach is not a real generalization of the genuine optimal transport, and by the way, to the best of our knowledge it was not implemented. Mass variations can manifest in different ways, we investigate first the case where two images have same maxima but still different mass. For instance consider the case where:

ρ 0 (x, y) = e -300 (x-0.3N ) 2 N 2 + (y-0.7N ) 2 2N 2 , ρ 1 (x, y) = e -200 (x-0.7N ) 2 2N 2 + (y-0.3N ) 2 N 2
where N × N is image size, we took N = 64 in the our tests. Note that the iso-contours are ellipses, of different orientations. We plot in figure 1 some pictures on the optimal path obtained with the exponential source term, and on figure 2(a) the mass variation during the corresponding interpolation. On figure 2(b) we depicted the

ρ(0) = ρ 0 ρ(0.25) ρ(0.5) ρ(0.75) ρ(1) = ρ 1 Fig 1.
Plot of the isolevels of the density ρ(t) along the optimal path between two densities of different mass.

residual decrease (in log scale) versus iterations, for the genuine no-source algorithm, the exponential and heat flux algorithms that we introduced. Note that the exponential algorithm performs very well, with a decreasing residual, on contrary with the no-source algorithm (that is beyond the domain of its validity, by the way). Another test consists in keeping the same support but consider different magnitudes. This is the case for the following initial and final densities, whose iso-contours are circles.

ρ 0 (x, y) = 1.5e As this is a symmetrical case, we depicted a side view on figure 3 to show the maximum decreasing, and in figure 4(a) the mass variation as well, still using our exponential source term. In figure 4(b) we plotted the residual decrease, and note that the exponential and heat flux source terms gave nearly the same decrease in that case. Note that in both cases, the optimal path is symmetric, that is, interpolating from ρ 0 to ρ 1 or from ρ 1 to ρ 0 gives the interpolating densities which are equal upon the transformation of t to 1t.

Next, we will consider the more intricate case where a mass splitting occurs. This could be relevant in tumor or more generally in cell growth. Our toy example corresponds to an initial gaussian function defined by As is clearly seen from figure 5, the genuine Benamou-Brenier algorithm attempts to create an optimal pas by adding mass at the very end of the path, while our algorithm produces a smooth mass variation which is more what we could expect form the interpolation of these two unbalanced densities.

Non-rigid registration for a growing lung metastasis

In this section some numerical tests on the non-rigid registration in biomedical imagery is proposed. Registration is adopted to refer to the same geometric configuration, especially when highly deformable organs are considered. The optimal transport is an objective way to provide a geometric transformation to this end. However, when tissue is evolving, as in tumor growth, the mass (measured considering the grey scale intensity) is not constant and a simple renormalization may provide unphysical mappings. A more accurate transformation could be obtained by considering continuous models of tumor growth. These models are often parametric and need a calibration, which may result in a costly process from a computational stand point. The proposed approach is a good tradeoff, allowing to get an interpolation between the images that is close to the real dynamics, for a computational cost that is lower if compared to that of inverse problems. Moreover, once the initial velocity potential has been obtained, an extrapolation can be done on a short time scale, that provides an approximated prognosis. Moreover in contrast with other models, existence and uniqueness of minimizers is proved. At last we point out that the proposed images are X-rays images, that is, the density of one pixel is directly related to absorption by X-rays, which in turn is higher in tumor cells. In that respect, images densities carry some biological informations and it make sense to use it to model biological phenomena.

The numerical experiments described hereafter concerns the non-rigid registration of portions of lung tissue with metastatic nodules. Several cases are considered, corresponding to different behaviors of the tumor growth. The first case considered is a sequence of 5 CT scan, shown in Fig. 6. They are a 2D slice of a 3D CT scan, represented in grey scale. The original images were segmented by using a threshold method, after having eliminated the bronchia and bronchiola structures: let the density be 0 ≤ ρ ≤ 1. If ρ ≤ 0.5 → ρ = 0, so that the image used for the computation has a density with a compact support (on the image, the black color corresponds to ρ = 0) . Remark 4. The Benamou-Brenier algorithm (and its alternatives proximal implementations), and therefore our algorithm, works for initial and final densities that may have compact supports (i.e. be identically zero on some part of the domain). Therefore we do not add any small quantity to raise densities in all the test cases proposed. The three model of source that we proposed were tested by taking the first and the last image of the sequence as initial and final densities respectively. In Fig. 7 the mass curve are compared to the real data: black circles are the masses of the segmented images, the black line is the mass of the interpolation obtained by an exponential source, the blue and the red one are the mass curves for the normal growth and the heat flux guided model. The real time scale for this evolution is 45 months: for the present work the time has been renormalized to [0, 1]; the same choice has been performed for all the testcases. The growth, in this case, is rather complex. Let us remark that the data are affected by noise and large errors may occur in the evaluation of the mass of the images after segmentation (even 30-40% of uncertainty). The three model proposed behaved differently and, despite the lack of biological modeling, the exponential and the normal sources are able to render the overall mass evolution with an error that is high but not larger than the measurement error, in some cases. The same considerations hold true for all the subsequent cases. In Fig. [START_REF] Benamou | Numerical Solution of the Optimal Transportation Problem via viscosity solutions for the Monge-Ampère equation[END_REF][START_REF] Dean | Numerical methods for fully nonlinear elliptic equations of the Monge-Ampère type[END_REF][START_REF] Figalli | A new transportation distance between non-negative measures, with applications to gradients flows with Dirichlet boundary conditions[END_REF] and the real data is shown,when the exponential model of growth is used. The images appears more regular in terms of shape with respect to the original one, but the accordance is good. In Fig. 9 the interpolated image is featured by a larger mass than the datum, and the error is at its maximum (as it can be checked by looking at Fig. 7).

Fast super-exponential growth

The second case we considered is a fast growth. A tissue portion is considered, in which a lung metastasis is growing. In Fig. 7.1.1 the sequence of three realistic images is shown. We tested three different model of sources, namely the exponential growth, the non-linear normal model and the heat-guided flux model. In particular, the optimal transport between the first and the last image was computed. In Fig. 11 the mass curve as function of time is shown for the three model sources compared to the data. There is a considerable increase of the tumor mass (about 6000%), that makes this case particularly challenging. For this case, the exponential source model was the more accurate one from a quantitative point of view (compare the three curves with the datum available at t = 0.5).

Logistic-type growth

The last case we considered is a tumor which initially grows in a rapid way and then undergoes a plateau type of evolution. The same analysis done for the previous cases is performed. In Fig. 12 the mass curve is shown for the different source models. The three models behave quite similarly in this case, and the more accurate one, that is, the one which is closer, in terms of mass, to the intermediate datum, is the normal growth one (blue curve in Fig. 12). None of the models is able to render a plateau-type of solution. A perspective might concern the setting up of a logistic model source to deal with this kind of growth. After having investigated the mass properties, the interpolation of the image is checked. In Fig. 13 are qualitatively well represented. The models proposed provide quite a good interpolation of tissues evolution, albeit their simplicity and the fact that they disregard the biology that governs the phenomena involved.

ρ(0) = ρ 0 ρ(0.5) ρ(1.0) = ρ 1 0 0.

Morphogenesis of kidneys

In this subsection a biological tree growth is considered. Several realistic images of the branching process occurring in the kidneys morphogenesis are available. As for the tumor growth case, there is no aim in investigating or understanding the biological and biophysical phenomena involved, but a realistic interpolation is sought. In Fig. 15 the mass curve for three different mappings is considered. The red line represents the mass curve when the interpolation between images (1 -3) is constructed. The blue and the black line are the mass curve for the mappings (1 -4) and (1 -5) respectively. From a quantitative standpoint only the mapping (1 -3) provides a good, realistic result. As said, the error is mainly due to the 2D nature of the images. The interpolation between the images (1 -5) is analyzed. The interpolated images (obtained by transport with the optimal mapping) are compared to the images of the sequence at corresponding time.

The qualitative agreement between the interpolated images and the original one is quite good (see Fig. [START_REF] Mccann | Polar factorization of maps on Riemannian manifolds[END_REF][START_REF] Peyre | Non-asymptotic equivalence between W 2 distance and H -1 norm[END_REF][START_REF] Piccoli | Generalized Wasserstein distance and its application to transport equations with source[END_REF] ). All the main geometrical features of the kidney tree are well reconstructed, even if only the first and the last image of the sequence are used to build the mapping.

Conclusion

In this paper we proposed a model of optimal transport between two unbalanced densities, which coincides with the usual transport in the balanced case. The model relies on the addition of a source term in the mass conservation equation. Three kinds of source terms are proposed and studied for application to tumor growth.

For two of them, we propose an existence result for the optimal transport problem, and give numerical evidence of its convergence and ability to deal with mass growth in various tests cases where the usual transport fails.

The key properties of our model are: a finite speed of motion of mass, and a the recovering of usual optimal transport for balanced densities, which was not the case of existing models for unbalanced densities [START_REF] Benamou | Numerical resolution of an "unbalanced" mass transfer problem[END_REF][START_REF] Piccoli | Generalized Wasserstein distance and its application to transport equations with source[END_REF][START_REF] Piccoli | A generalized Benamou-Brenier formula for mass-varying densities[END_REF].

Proposition 4 .

 4 Let us consider the following density pair, for a smooth and positive function ρ 0 :

  c[ρ] as above and find c[ρ] = c, thus (ρ, m) ∈ C 2 , which is nonempty. Now take (ρ, m) ∈ C 2 , c[ρ] is a constant and integrating the ODE (60) give c[ρ] = c, thus (ρ, m) ∈ C 1 .

Fig 2 .= ρ 1 Fig 3 .

 213 Fig 2. (a) Plot of the mass variation from initial to final density during the optimal path for the exponential source term. (b) Plot of the residual variation during iterations, for the no-source, exponential source, and heat flux source algorithms. Note that the vertical axis is log-scaled.

Fig 4 .

 4 Fig 4. (a) Plot of the mass variation from initial to final density during the optimal path of the exponential source term. (b) Plot of the residual variation during iterations, for the no-source, exponential source, and heat flux source algorithms. Note that the vertical axis is log-scaled.

2 += ρ 1 Fig 5 .

 215 Fig 5. Plot of the isolevels of the density ρ(t) along the optimal path between two densities of different mass. Top rows: original Benamou-Brenier Algorithm, with top and side views; Bottom rows: Our algorithm with exponential source term.

= ρ 1 ρ 1 Fig 6 .

 116 Fig 6. Plot of a growing tumor simulation from an initial and final stage. Top: B/W images; Bottom: corresponding density plots.

Fig 7 .

 7 Fig 7. Mass curve as function of time for different source models (color continuous line), compared to the data (black circles).

Fig 8 .Fig 9 .

 89 Fig 8. Comparison between the interpolation at time t = 0.5, a), obtained by solving optimal transport between images (1 -5) and the original image in grey scale b)

Fig 10 .Fig 11 .

 1011 Fig 10. Comparison between the interpolation at time t = 0.875, a), obtained by solving optimal transport between images (1 -5) and the original image in grey scale b)

Fig 12 .

 12 Fig 12. Mass curve as function of time for different source models (color continuous line), compared to the data (black circles).

Fig 13 .

 13 Fig 13. Comparison between the interpolation at time t = 0.8, a), obtained by solving optimal transport between images (1 -3) and the original image in grey scale b)

1 Fig 14 .Fig 15 .

 11415 Fig 14. Images of kidneys morphogenesis

Fig 16 .Fig 17 .Fig 18 .

 161718 Fig 16. Comparison between the interpolation at time t = 0.25, a), obtained by solving optimal transport between images (1 -5) and the original image in grey scale b)

imsart-generic ver. 2012/08/31 file: dlemOrth.tex date: October 7, 2014

Acknowledgment: We would like to thank Q. Mérigot who pointed out to us the references [18,19]. E. Maitre is supported by French Agence Nationale de la Recherche, ANR Project TOMMI (ANR 2011 BS01 014 01) and by Grenoble INP, through SEI grant TOSCANA.