Nicolas Papadakis

Etienne Mémin

Juin

A variational method for joint tracking of curve and motion

Keywords: Optimal control, dynamical model, tracking of curves, variational method

In this paper, a new framework for the tracking of closed curves and their associated motion fields is described. The proposed approach enables a continuous tracking along an image sequence of both a deformable curve and velocity field. Such an approach is formalized through the minimization of a global spatio-temporal continuous cost functional, w.r.t a set of variables representing the curve and its related motion field. Relying on an optimal control technique, the resulting minimization sequence consists in a forward integration of an evolution law followed by a backward integration of an adjoint evolution model. This latter pde includes a term related to the discrepancy between the state variables evolution law and discrete noisy measurements of the system. The closed curves are represented through implicit surface modeling, whereas the motion is described either by a vector field or through vorticity and divergence maps depending on the kind of targeted applications.

The efficiency of the approach is demonstrated on two types of image sequences showing deformable objects and fluid motions.

Mots-clés : Contrôle optimal; modèle dynamique; suivi de courbes; méthode variationnelle 1 Introduction

Motivation and scope

Tracking the contours and the motion of an object is an essential task in many applications of computer vision. For several reasons such a generic issue appears to be very challenging in the general case. As a matter of fact, the shape of a deformable object or even of a rigid body may change drastically when visualized from an image sequence. These deformations are due to the object's proper apparent motion or to perspective effect and 3D shape evolution.

This difficulty is amplified when the object becomes partially or totally occluded during even a very short time period. Such a visual tracking is also a challenging issue when the curve of interest delineates iso-quantities transported by a fluid motion. This last case is of importance in domains such as meteorology or oceanography where one may wish to track iso-temperature, contours of cloud systems, or the vorticity of a motion field. Here, the most difficult technical aspect consists in handling the tracking of these features in a consistent way with respect to appropriate physical conservation laws. Another serious difficulty comes from the dimensionality of the state variables. In contexts involving for instance geophysical flows, the targeted curve which represents clouds, sea color or temperature may exhibit high topological changes and presents unsteady irregularities over time. Such behaviors make difficult the use of a reduced parametric description of the curves. They can be only described in spaces of very large dimension. In addition, the fluid motion of the regions enclosed by such curves may only be accurately described by dense motion fields. The joint tracking of a curve and its underlying motion requires as a result to handle a state space of huge dimension. This curse of dimensionality remains true for general unknown natural objects (either rigid or deformable) observed in complex environments.

Related works and their limitations

This context makes difficult the use of recursive Bayesian filters such as the particle filter [START_REF] Blake | Active contours[END_REF],

since stochastic sampling in large state spaces is usually completely inefficient. Even if such a filter has recently been used for fluid motion tracking [START_REF] Cuzol | A low dimensional fluid motion estimator[END_REF][START_REF] Cuzol | A stochastic filter for fluid motion tracking[END_REF], this kind of techniques are only applicable when the unknown state can be described on a reduced set of basis functions. Furthermore, coping with a coupled tracking of curves and motion fields augments significantly the problem complexity. For such an issue, numerous approaches based on the level set representation have been proposed [START_REF] Cremers | Motion competition: A variational framework for piecewise parametric motion segmentation[END_REF][START_REF] Goldenberg | Fast geodesic active contours[END_REF][START_REF] Kimmel | Tracking level sets by level sets: a method for solving the shape from shading problem[END_REF][START_REF] Mansouri | Region tracking via level set PDEs without motion computation[END_REF][START_REF] Niethammer | Dynamic geodesic snakes for visual tracking[END_REF][START_REF] Osher | Fronts propagating with curvature dependent speed: Algorithms based on Hamilton-Jacobi formulation[END_REF][START_REF] Paragios | Geodesic active regions: a new framework to deal with frame partition problems in computer vision[END_REF][START_REF] Sethian | Level set methods[END_REF]. All these techniques describe the tracking as successive 2D segmentation processes sometimes enriched with a motion based propagation step. Segmentation techniques on spatio-temporal data have also been proposed [START_REF] Bartesaghi | Tracking of moving objects under severe and total occlusions[END_REF][START_REF] Goldenberg | Fast geodesic active contours[END_REF]. Since level sets methods do not introduce any temporal consistency related to a given dynamical laws -i.e. a tracking process -, they are quite sensitive to noise [START_REF] Martin | Influence of the noise model on level set active contours segmentation[END_REF] and exhibit inherent temporal instabilities. Implausible growing/decreasing or merging/splitting cannot

Papadakis & Mémin

be avoid without introducing hard adhoc constraints or some statistical knowledges on the shape [START_REF] Cremers | Diffusion snakes: introducing statistical shape knowledge into the Mumford-Shah functional[END_REF][START_REF] Leventon | Statistical shape influence in geodesic active contours[END_REF][START_REF] Paragios | A level set approach for shape-driven segmentation and tracking of the left ventricle[END_REF] and as a consequence dedicate the process to very specific studies. Besides these approaches can hardly handle occlusions of the target or to cope with severe failures of the image measurements (for instance a complete loss of image data, a severe motion blur, high saturation caused by over exposure or failure of the low level image detectors). In the context of fluid motion, a prior learning of the "fluid object" shape is in essence almost impossible and the inclusion of fluid dynamical laws related to the Navier-Stokes equation is essential to provide coherent physical plausible solution.

In [START_REF] Soatto | Geometric level sets methods in imaging vision and graphics[END_REF], an approach based on a group action mean shape has been used in a moving average context. Contrary to previous methods, this approach introduces, through the moving average technique, a kind of tracking process. This tracking is restricted to simple motions and does not allow to introduce complex dynamical law defined through differential operators.

The explicit introduction of a dynamic law in the curve evolution law has been considered in [START_REF] Niethammer | Dynamic geodesic snakes for visual tracking[END_REF]. However, the proposed technique needs a complex detection mechanism to cope with occlusions.

Contribution

In this paper, we propose a technique which allows to handle both the tracking of closed curves and the underlying motion field transporting this curve. The approach is related to variational data assimilation technique used for instance in meteorology [START_REF] Bennet | Inverse Methods in Physical Oceanography[END_REF][START_REF] Dimet | Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects[END_REF][START_REF] Talagrand | Variational assimilation of meteorological observations with the adjoint vorticity equation[END_REF]. Such a technique enables, in the same spirit as a Kalman filter, a temporal smoothing along the whole image sequence. It combines a dynamical evolution law of state variables representing the target of interest with the whole set of available measurements related to this target.

Unlike Bayesian filtering approach which aim at estimating a probability law of the stochastic process associated to the feature of interest, variational assimilation has the advantage to allow handling state spaces of high dimension. From the motion analysis point of view, such framework allows to incorporate a dynamical consistency along the image sequence.

The framework we propose provides an efficient technique to incorporate general dynamical constraints into a coupled motion and segmentation process. As it is demonstrated in this paper, the method is particularly well suited to the analysis of fluid motion and deformations. This technique could also be useful for batch video processing. We provide some results toward this direction even if efficient answers for such applications would need much more works.

Outline of the paper

The paper is organized as follows. After a description of the data assimilation technique in section 2, we introduce the proposed curve tracking method in section 3. The method is then extended for a joint motion and object tracking in section 4. Two kind of applications are studied here: fluid motion tracking in section 4.1 and natural object tracking in section 4.2. As it will be demonstrated in the experimental sections, such a technique enables to handle naturally noisy data and complete loss of image data over long time periods without resorting to complex mechanisms. This paper extends a conference paper which focus on contours tracking [START_REF] Papadakis | Variational optimal control technique for the tracking of deformable objects[END_REF].

Variational tracking formulation

In this section, we first describe the general framework proposed for tracking problems.

It relies on variational data assimilation concepts [START_REF] Bennet | Inverse Methods in Physical Oceanography[END_REF][START_REF] Dimet | Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects[END_REF][START_REF] Talagrand | Variational assimilation of meteorological observations with the adjoint vorticity equation[END_REF] proposed for the analysis of geophysical flows. For a sake of clarity, we will first present a simplified formulation where the state variable of interest obeys to a perfect evolution law. Let us note that this situation corresponds to the most usual model used in geophysical domain for data assimilation. This section will also be for us a mean to define properly the basic ingredients of the framework.

Data Assimilation with perfect model

Direct evolution model Let the state space V be an Hilbert space identified to its dual space. Noting X ∈ W(t 0 , t f ) the state variable representing the feature of interest, which is assumed to live in a functional space

W(t 0 , t f ) = {X|X ∈ L 2 (t 0 , t f ; V), ∂ t X ∈ L 2 (t 0 , t f ; V)}
and assuming that the evolution in time range [t 0 ; t f ] of the state is described through a (non linear) differential model M : V×]t 0 , t f [→ V, we get the following direct problem:

For a given η ∈ V, let us define X ∈ W(t 0 , t f ) such that:

∂ t X(t) + M(X(t), t) = 0, X(t 0 ) = X 0 + η. (1)
This system gathers an evolution law and an initial condition of the state variable. It is governed by a control variable η ∈ V, identified here to the inaccuracy on the initial condition.

The control could also be defined on model's parameters [START_REF] Dimet | Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects[END_REF]. The direct problem (1) will be assumed to be well posed, which means that we first assume that the application V → V : η → X(η, t) is differentiable ∀t ∈]t 0 , t f [ and secondly that given η ∈ V and t f > t 0 , there exists a unique function X ∈ W(t 0 , t f ) solution of problem [START_REF] Anderson | Optimal Filtering[END_REF] and that this solution depends continuously on η (i.e:

V → V : η → X(η, t) is continuous ∀t ∈]t 0 , t f [).
Let us also assume that some measurements (also called observations) Y ∈ O of the state variable components are available. These observations may live in a different space (a reduced space for instance) from the state variable. We will nevertheless assume that there exists a (non linear) observation operator H : V → O, that goes from the variable space to the observation space.

Papadakis & Mémin

Cost function Let us define an objective function J : V → R measuring the discrepancy between a solution associated to an initial state control variable of high dimension and the whole sequence of available observations as:

J(η) = 1 2 t f t0 ||Y -H(X(η, t), t)|| 2 R dt + 1 2 ||η|| 2 B . (2)
The overall problem that we are facing consists in finding the control variable η ∈ V that minimizes the cost function J. Norms || • || R and || • || B are respectively associated to the

scalar products R -1 •, • O and B -1 •, • V ,
where R and B are symmetric positive defined endomorphisms of V. In our applications, R and B are respectively called, with some abuse of language, the observation covariance matrix and the initialization covariance matrix. They enable to weight the deviations from the true initial state to the given initial condition and to quantify the exactitude of the observation model.

Differential model

In order to compute the partial derivative of the cost function with respect to the control variable, system (1) is differentiated with respect to η in the direction δη. The following differential model is obtained:

Given η ∈ V , X(t) a solution of (1) and a perturbation δη ∈ V,

dX = ∂X ∂η δη ∈ W(t 0 , t f ) is such that: ∂ t dX(t) + (∂ X M)dX(t) = 0, dX(t 0 ) = δη. (3)
In this expression, the tangent linear operator (∂ X M) is defined as the Gâteaux derivative of the operator M at point X:

(∂ X M)dX(t) = lim β→0 M(X(t) + βdX(t)) -M(X(t)) β . (4)
The tangent linear operator (∂ X H) associated to H may be defined similarly. Differentiating now the cost function (2) with respect to η in the direction δη leads to:

∂J ∂η , δη V = - t f t0 Y -H(X), (∂ X H) ∂X ∂η δη R dt + X(t 0 ) -X 0 , δη B = - t f t0 (∂ X H) * R -1 (Y -H(X)), ∂X ∂η δη V dt + B -1 (X(t 0 ) -X 0 ), δη V . (5)
where (∂ X H) * , the adjoint operator of (∂ X H), is defined by the scalar product:

∀x ∈ V, ∀y ∈ O (∂ X H) x, y O = x, (∂ X H) * y V . (6)
Adjoint evolution model In order to estimate the gradient of the cost function J, a first numerical brute force approach consists in computing the functional gradient through finite differences:

∇ u J ≃ J(η + ǫe k ) -J(η) ǫ , k = 1, . . . , p ,
where ǫ ∈ R is an infinitesimal perturbation and {e k , k = 1, . . . , p} denotes the unitary basis vectors of the control space V. Such a computation is impractical for space of large dimension since it requires p integrations of the evolution model for each required value of the gradient functional. Adjoint models as introduced in optimal control theory by J.L. Lions [START_REF] Lions | Contrôle optimal de systèmes gouvernés par des équations aux dérivées partielles[END_REF][START_REF] Lions | Optimal control of systems governed by partial differential equations[END_REF] and seminally applied in meteorology in [START_REF] Dimet | Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects[END_REF] will allow us to compute the gradient functional in a single integration. To obtain the adjoint equation, the first equation of model ( 3) is multiplied by an adjoint variable λ ∈ W(t 0 , t f ). The result is then integrated on [t 0 , t f ]:

t f t0 ∂ t dX(t), λ(t) V dt + t f t0 (∂ X M)dX(t), λ(t) V dt = 0,
After an integration by parts of the first term and using the second equation of the differential model (3), we finally get:

-

t f t0 -∂ t λ(t) + (∂ X M) * λ(t), dX(t) V dt = λ(t f ), dX(t f ) V -λ(t 0 ), δη V , (7)
where the adjoint operator (∂ X M) * is defined by the scalar product:

∀x ∈ V, ∀y ∈ V (∂ X M) x, y V = x, (∂ X M) * y V . (8)
In order to obtain an accessible expression for the cost function gradient, we impose that λ(t f ) = 0 and we define the following adjoint problem:

Given η ∈ V, t f > t 0 and X(t) solution of (1), let us define λ ∈ W(t 0 , t f ) such that:

-∂ t λ(t) + (∂ X M) * λ(t) = (∂ X H) * R -1 (Y -H(X(t))) ∀t ∈]t 0 , t f [, λ(t f ) = 0. (9)
In the same way as for the direct model, we assume that given η ∈ V, t f > t 0 and X ∈ W(t 0 , t f ) solution of problem (1), there exists a unique function λ ∈ W(t 0 , t f ) solution of problem [START_REF] Chan | Active contours whitout edges[END_REF]. We also assume that this solution depends continuously on η (i.e:

V → V : η → λ(η, t) is continuous ∀t ∈]t 0 , t f [).
Functional gradient Combining now equations ( 5), ( 7) and ( 9), we obtain:

∂J ∂η , δη V = -λ(t 0 ), δη V + B -1 (X(t 0 ) -X 0 ), δη V .
The cost function derivative with respect to η finally reads:

∂J ∂η = -λ(t 0 ) + B -1 (X(t 0 ) -X 0 ). ( 10
)
As a consequence, given a solution X(t) of the direct model (1), the functional gradient can be computed with a backward integration of the adjoint model [START_REF] Chan | Active contours whitout edges[END_REF]. The adjoint variable then enables to update the initial condition, by canceling the gradient defined in [START_REF] Corpetti | Dense estimation of fluid flows[END_REF]:

Papadakis & Mémin X(t 0 ) = X 0 + Bλ(t 0 ) ( 11
)
where B is the pseudo inverse of B -1 [START_REF] Bennet | Inverse Methods in Physical Oceanography[END_REF]. As the complexity of the adjoint model integration is similar to the integration of the direct model, the use of this technique appears to be very efficient for state space of large dimension. A synoptic of the overall technique is given in algorithm (2.1). This first approach is widely used in environmental sciences for the Algorithm 2.1 Perfect Model

Let X(t0) = X0.

(i) From X(t0), compute X(t), ∀t ∈]t0, t f [ with a forward integration of system (1).

(ii) With X(t), realize a backward integration of the adjoint variable with the system (9).

(iii) Update the initial condition X(t0) with relation [START_REF] Courtier | Variational assimilation of meteorological observations with the direct and adjoint shallow-water equations[END_REF].

(iv) Return to (i) and repeat until a convergence criterion.

analysis of geophysical flows [START_REF] Dimet | Variational algorithms for analysis and assimilation of meteorological observations: theoretical aspects[END_REF][START_REF] Talagrand | Variational assimilation of meteorological observations with the adjoint vorticity equation[END_REF]. For these applications, the involved dynamical model are assumed to reflect faithfully the evolution of the observed phenomenon. However, such modelling seems to us irrelevant in image analysis since the different models on which we can rely on are usually inaccurate due for instance to 3D-2D projections, varying lighting conditions, or completely unknown boundary conditions. Considering imperfect dynamical models now comes to an optimization problem where the control variable is related to the whole trajectory of the state variable. This is the kind of problem we are facing in this work.

Data Assimilation with imperfect model

The dynamical model we consider now is defined up to a control function ν ∈ W(t 0 , t f , V),

where ν(t) ∈ V, and as previously up to a control variable on the initial state η ∈ V. We are now facing an imperfect dynamical system which depends on the whole trajectory of the model control function and on the initial state control variable. Formally, the system associated to an imperfect model reads:

Given (ν, η) ∈ (W, V), let us define X ∈ W(t 0 , t f ) such that ∂ t X(t) + M(X(t), t) = ν(t) ∀t ∈]t 0 , t f [, X(t 0 ) = X 0 + η. ( 12
)
Cost function As previously, the objective function J : W ×V → R gathers a measurement discrepancy term and penalization terms on the control variables norms:

J(ν, η) = 1 2 t f t0 ||Y -H(X(ν(t), η, t))|| 2 R dt + 1 2 ||η|| 2 B + 1 2 t f t0 ||ν(t)|| 2 Q dt. (13)
The norm || • || Q is associated to the scalar product Q -1 •, • V , where Q is a symmetric positive defined endomorphism of V called the model covariance matrix. We aim here at finding a minimizer (η, ν) of the cost function that is of least energy and that minimizes along time the deviations between the available measurements and the state variable.

Differential model In the aim of computing partial derivative of the cost function with respect to the control variables, system ( 12) is first differentiated with respect to (ν, η) in the direction (δν, δη):

Given (ν, η)∈(W, V), X(t) solution of ( 12) and a perturbation (δν, δη)∈(W × V),

dX = ∂X ∂ν δν + ∂X ∂η δη ∈ W(t 0 , t f ), is such that: ∂ t dX(t) + (∂ X M)dX(t) = δν(t) ∀t ∈]t 0 , t f [, dX(t 0 ) = δη. (14)
The differentiation of cost function [START_REF] Cremers | Diffusion snakes: introducing statistical shape knowledge into the Mumford-Shah functional[END_REF] with respect to η has been previously computed and given in equation ( 5). The differentiation with respect to ν in the direction δν reads:

∂J ∂ν , δν W = - t f t 0 Y -H(X), (∂X H) ∂X ∂ν δν(t) R dt + t f t 0 ∂tX(t) + M(X(t)), δν(t) Q dt = - t f t 0 (∂XH) * R -1 (Y -H(X)), ∂X ∂ν δν(t) V dt + t f t 0 Q -1 (∂tX(t)+M(X(t)), δν(t) V dt. (15)
Adjoint evolution model Similarly to the previous case, the first equation of model ( 14) is multiplied by an adjoint variable λ ∈ W(t 0 , t f ) and integrated on [t 0 , t f ]

t f t 0 ∂tdX(t), λ(t) V dt + t f t 0 ∂XM(X(t))dX(t), λ(t) V dt = t f t 0 δν(t), λ(t) V dt.
After an integration by parts of the first term and using the second equation of the differential model ( 14), we finally get:

-

t f t 0 -∂tλ(t) + (∂XM) * λ(t), dX(t) V dt = λ(t f ), dX(t f ) V -λ(t0), δη V - t f t 0 λ(t), δν(t) V dt. (16)
As previously, in order to exhibit an expression of the gradient of the cost function from the adjoint variable, we define the following adjoint problem:

Given (ν, η) ∈ (W, V), t f > t 0 and X(t) a solution of (12), we define λ ∈ W(t 0 , t f ), such that:

-∂ t λ(t) + (∂ X M) * λ(t) = (∂ X H) * R -1 (Y -H(X(t))) ∀t ∈]t 0 , t f [, λ(t f ) = 0. ( 17
)
Functional gradient Combining equations ( 5) and [START_REF] Cuzol | A stochastic filter for fluid motion tracking[END_REF], the functional gradient is given by:

∂J ∂ν , δν W + ∂J ∂η , δη V = t f t 0 Q -1 (∂tX(t)+M(X(t)), δν(t) V dt + B -1 (X(t0)-X0), δη V - t f t 0 (∂X H) * R -1 (Y -H(X(t))), ∂X ∂w δν(t) + ∂X ∂η δη dX(t) p V dt.
Introducing ( 16) and ( 17), we obtain:

∂J ∂ν , δν W + ∂J ∂η , δη V = t f t 0 Q -1 (∂tX(t)+M(X(t))-λ(t), δν(t) V dt-λ(t0), δη V + B -1 (X(t0) -X0), δη V = Q -1 (∂tX + M(X) -λ, δν W + -λ(t0) + B -1 (X(t0) -X0), δη V .
The derivatives of the cost function with respect to ν and η are identified as:

∂J ∂ν = Q -1 (∂tX + M(X)) -λ, (18)
∂J ∂η

= -λ(t0) + B -1 (X(t0) -X0). (19)
Canceling these components and introducing Q and B, the respective pseudo inverses of Q -1 and B -1 , we get:

∂tX(t) + M(X(t)) = Qλ(t) X(t0) -X0 = Bλ(t0). ( 20
)
The second equation still constitutes an incremental update of the initial condition. However, the integration of this system requires the knowledge of the whole adjoint variable trajectory, which itself depends on the state variable through [START_REF] Faugeras | Assimilation variationnelle de donnes dans un modle coupl ocanbiogochimie[END_REF]. This system can be in practice integrated introducing an incremental splitting strategy.

Incremental function Defining the state variable as

X(t) = X(t) + dX(t) ∀t ∈ [t0, t f ], X(t0) = X0, (21)
where X(t) is either a fixed component or a previously estimated trajectory of the state variable, equation ( 20) can be written as:

∂t X(t) + M( X(t)) = 0 ∀t ∈]t0, t f [, (22)
∂tdX(t) + ∂ X M( X(t))dX(t) = Qλ(t) ∀t ∈]t0, t f [. ( 23
)
The update of the state variable X is as a consequence driven by an incremental function dX.

The adjoint variable λ is obtained from a backward integration of ( 17) and the trajectory of X. The initial value of this incremental function is otherwise obtained from [START_REF] Harten | Uniformly high order accurate essentially non-oscillatory schemes, 111[END_REF] and reads:

dX(t 0 ) = Bλ(t 0 ). (24)
Equations ( 12), ( 21), ( 22), ( 23) and ( 24) give rise to a variational assimilation method with imperfect dynamical model. A sketch of the whole process is summarized in Algorithm (2.2).

The algorithm principles are also schematically pictured on figure 1.

Algorithm 2.2 Let X(t 0 ) = X 0 .

(i) From X(t 0 ), compute X(t), ∀t ∈]t 0 , t f [ with a forward integration of system [START_REF] Kohlberger | Variational dense motion estimation using the helmholtz decomposition[END_REF].

(ii) X(t) being given, realize a backward integration of the adjoint variable with the system [START_REF] Faugeras | Assimilation variationnelle de donnes dans un modle coupl ocanbiogochimie[END_REF].

(iii) Compute the initial value of the incremental function [START_REF] Kurganov | New high-resolution central schemes for nonlinear conservation laws and convection-diffusion equations[END_REF].

(iv) From dX(t 0 ), compute dX(t), ∀t ∈]t 0 , t f [ with a forward integration of system ( 23).

(v) Update X = X + dX.

(vi) Return to (ii) and repeat until convergence. of the adjoint variable relying on a measurement discrepancy (purple trajectory) enables to compute a forward incremental correction trajectory (red curve) and so on...

Additional ingredients

Before turning to the application of such a framework, let us note that the system we have presented so far can be slightly modified and enriched considering a final condition or additional types of observations.

Final condition Symmetrically to the initial condition equation, a final target state can eventually be added through an additional equation:

Given (ν, η, w) ∈ (W, V, V), let us define X ∈ W(t 0 , t f ), such that:        ∂ t X(t) + M(X(t), t) = ν(t) ∀t ∈]t 0 , t f [, X(t 0 ) = X 0 + η, X(t f ) = X f + w, (25)
with w ∈ V an additional control variable. The cost function to minimize in this case incorporates an extra penalization term on the norm of this new control variable:

J(ν, η, w) = 1 2 t f t 0 ||Y -H(X(ν(t), η, w, t))|| 2 R dt + 1 2 ||η|| 2 B + 1 2 t f t 0 ||ν(t)|| 2 Q dt + 1 2 ||w|| 2 F . (26)
The norm || • || F is associated to the scalar product F -1 •, • I . The differentiation of the cost function [START_REF] Leventon | Statistical shape influence in geodesic active contours[END_REF] with respect to ν in the direction δν is:

∂J ∂w , δν V = - t f t0 (∂ X H) * R -1 (Y -H(X)), ∂X ∂w δw V dt + F -1 (X(t f ) -X f ), δw V . (27)
The introduction of an additional equation on the final state modifies the final condition of the adjoint model. It now reads:

Given (ν, η, w) ∈ (W, V, V), t f > t 0 and X(t) a solution of ( 25), let the adjoint variable λ ∈ W(t 0 , t f ), be defined such that:

-∂ t λ(t) + (∂ X M) * λ(t) = (∂ X H) * R -1 (Y -H(X)) ∀t ∈]t 0 , t f [, λ(t f ) = F -1 (X f -X(t f )). ( 28
)
Several observations If we consider a set of observations Y i ∈ O i , i ∈ {1, . . . , N } related to the state variable through the observation operators H i (V, O i ), the measurement discrepancy term of the cost function becomes:

1 2 N i=1 ||Y i -H i (X, η, t)|| 2 R i . (29)
The norms || • || R i involved in this term are associated to the scalar products

R i -1 •, • Oi .
Such a term also implies a straightforward modification of the adjoint model:

Given η ∈ V, t f > t 0 and X(t) solution of (1), let the adjoint λ ∈ W(t 0 , t f ), be such that:

-∂ t λ(t) + (∂ X M) * λ(t) = N i=1 (∂ X H i ) * R -1 i (Y i -H i (X)) ∀t ∈]t 0 , t f [ λ(t f ) = 0. ( 30
)

Discussion on convergence

The convergence criterion introduced in Algorithm 2.1 can be defined from the norm of the objective function gradient, namely |λ(t 0 )| < α, where α is a given threshold. The value of this threshold must be as small as possible, but above numerical errors due to the discretization schemes used. An empirical test related to the numerical stability of the adjoint operators discretization [START_REF] Courtier | Variational assimilation of meteorological observations with the direct and adjoint shallow-water equations[END_REF][START_REF] Faugeras | Assimilation variationnelle de donnes dans un modle coupl ocanbiogochimie[END_REF] allows to determine this threshold value. Referring to Taylor expansion, this test consists in checking that the ratio

lim α→0 J(X + αdX) -J(X) α∇J → 1.
Computing the numerator through finite differences, whereas the denominator is obtained by the adjoint model, the test enables to find for a given set of measurements the smallest value of α that respects this limit.

In order to optimize the minimization process, different types of efficient gradient descent strategies can be used (Conjugate Gradient, quasi Newton methods...). In this work, we use a simple fixed step gradient descent.

As for the global convergence of the method, the minimization problem has a unique solution if the functional J is convex, lower semi-continuous and if:

lim ||ν||W →∞,||η||V →∞ J(ν, η) = ∞.
If the involved models M and H are non linear, the functional will be likely not convex. The convergence (assuming that there exists a unique solution which depends continuously on the control variable) may be realized only toward a local minimum which depends on the chosen initialization. Hence, the initialization term is, in general, of great importance to obtain a good solution. However, the control parameter on the initial condition equation enables to model an incertitude on this initial state.

Relations with the Kalman filter and the Kalman smoother

We briefly discuss in this section the relations between the previous system and the Kalman filter and smoother. As variational assimilation technique and Kalman filtering are based on very same kind of systems, both techniques share some similarities. Theoretical equivalence between the two techniques have been proved for linear systems (both model and observation operators are linear) [START_REF] Bouttier | Data assimilation concepts and methods[END_REF][START_REF] Li | Optimality of 4D-var and its relationship with the kalman filter and kalman smoother[END_REF]. One of the main difference concerns the goal sought by each technique. Kalman filter aims at computing recursively in time the two first moments of the conditional distribution of the state variables given all the available observations [START_REF] Anderson | Optimal Filtering[END_REF] formally the sequence of past observations in case of filtering and the complete sequence in case of smoothing -whereas a variational assimilation technique aims at estimating the minimal cost trajectory in batch mode. In case of a linear Gaussian system, the first moment of the filtering distribution and the solution of the variational technique are equivalent. This is not the case for non linear systems. Furthermore, unlike Kalman filtering, variational approaches do not provide any covariance of the estimation error. From a practical and computational point of view, there exists also some differences between both techniques.

Kalman filtering requires at each iteration to inverse the estimation error covariance matrix in order to compute the so called Kalman gain matrix. The dimension of this matrix is the square of the size of the state vector, as a consequence Kalman filtering is unfeasible for large state vectors. At the opposite, variational assimilation does not require such an inversion. The estimation is done here iteratively by forward-backward integrations and not directly through a recursive process like in Kalman filtering. This point makes variational assimilation methods very attractive for the tracking of high dimensional features.

Batch filtering vs. recursive tracking

As already mentioned in the previous section, one of the main differences between Bayesian filters and variational assimilation techniques relies in the underlying integration that is operated. The forward recursive expression of the former makes Bayesian filter well suited to real time tracking. They are on the other side only efficient for low dimensional state spaces.

The latter techniques are inherently batch processes as they are based on forward-backward integration schemes. Unless relying on temporal sliding windows such a characteristic makes their use difficult for real time tracking. Nevertheless, their abilities to cope with large dimensional state spaces and with (highly) non linear differential dynamics make them interesting for a batch video processing. Such analysis is interesting each time that one wishes to extract from a sequence of noisy and incomplete data a continuous and temporal coherent sequence of features with respect to a specified dynamics. As we will show in the following sections this is of particular interest when one aims at analyzing fluid flows from image sequences. We also believe that such a scheme could bring valuable information for the analysis of deformations in medical images or in the domain of video processing for purposes of video inpainting, object recolorization, morphing or video restoration. We give some hints about these potential applications in the following.

Application to curve tracking

Before dealing with the most complex issue of jointly tracking curve and motion, we will first focus in this section on the application of the variational data assimilation to curve tracking from an image sequence. In this first application, the motion fields driving the curve will be considered as reliable external inputs of the system. For such a system, different types of measurements will be explored.

Contour representation and evolution law

As we wish to focus in this work on the tracking of non parametric closed curves that may exhibit topology changes during the time of the analyzed image sequence, we will rely on an implicit level set representation of the curve of interest Γ(t) at time t ∈ [t 0 , t f ] of the image sequence [START_REF] Osher | Fronts propagating with curvature dependent speed: Algorithms based on Hamilton-Jacobi formulation[END_REF][START_REF] Sethian | Level set methods[END_REF]. Within that framework, the curve Γ(t) enclosing the target to track is implicitly described by the zero level set of a function φ(x, t) : Ω × R+ → R:

Γ(t) = {x ∈ Ω | φ(x, t) = 0},
where Ω stands for the image spatial domain. This representation enables an Eulerian representation of the evolving contours. As such, it allows to avoid the ad-hoc regridding processes of the different control points required in any explicit Lagrangian -spline based -curve description. The problem we want to face consists in estimating for a whole time range the state of an unknown curve, and hence of its associated implicit surface φ. To that end, we first define an a priori evolution law of the unknown surface. We will assume that the curve is transported at each frame instant by a velocity field, w(x, t) = [u(x, t), v(x, t)] T , and diffuses according to a mean curvature motion. We will assume that this evolution law is verified only up to an additive control function. In term of implicit surface, the overall evolution law reads:

∂ t φ + ∇φ • w -εκ ∇φ = ν, (31)
where κ denotes the curve curvature, and ν(x, t) the control function. Introducing the surface normal, equation ( 31) can be written as:

∂ t φ + (w • n -εκ) ∇φ △ =M(φ) = ν, (32)
where the curvature and the normal are directly given in term of surface gradient:

κ = div ∇φ ∇φ and n = ∇φ ∇φ .
As previously indicated, the motion field transporting the curve will be first assumed to be given by an external estimator. In practice, we used an efficient and robust version of the

Papadakis & Mémin

Horn and Schunck optical-flow estimator [START_REF] Mémin | Dense estimation and object-based segmentation of the optical flow with robust techniques[END_REF]. Concerning the covariance Q of the control term involved in the evolution law, we fixed it to a constant diagonal matrix since it is rather difficult to infer precise model's errors. Two different ranges of value have been used. For sequences where accurate motion fields can be recovered we fixed this value to a low value (typically 0.005). At the opposite, for sequences for which only motion field of bad quality are obtained this parameter has been set to a higher value (0.5) as we are facing in that case a rough dynamic.

Tangent linear evolution operator To apply the variational data assimilation principles, we must first define the expression of the Gâteaux derivative of the operators involved. Using equation ( 32), the evolution operator reads in its complete form as:

M(φ) = ∇φ • w -ε||∇φ||div ∇φ ||∇φ|| .
This operator can be turned into a more tractable expression for our purpose:

M(φ) = ∇φ T w -ε ∆φ - ∇ T φ∇ 2 φ∇φ ||∇φ|| 2 . ( 33
)
The corresponding tangent linear operator at point φ finally reads:

∂ φ Mδφ = ∇δφ T w -ε ∆δφ - ∇φ T ∇ 2 δφ∇φ ||∇φ|| 2 + 2 ∇φ T ∇ 2 φ ||∇φ|| 2 ∇φ∇φ T ||∇φ|| 2 -Id ∇δφ . (34)
Operator discretization Before going any further, let us describe the discretization schemes we considered for the evolution law. This concerns the evolution operator, the associated tangent linear operator and the adjoint evolution operator. We will denote by φ t i,j the value of φ at image grid point (i, j) at time t ∈ [t 0 ; t f ]. Using [START_REF] Lions | Optimal control of systems governed by partial differential equations[END_REF] and a semi-implicit discretization scheme, the following discrete evolution model is obtained:

φ t+∆t i,j -φ t i,j ∆t + M φ t i,j φ t+∆t i,j = 0.
Considering φ x and φ y , the horizontal and vertical gradient matrices of φ, the discrete operator M is obtained as :

M φ t i,j φ t+∆t i,j = (φ t+∆t x )i,j (φ t+∆t y )i,j T w - ε ||∇φ t i,j || 2 -(φ t y )i,j (φ t x )i,j T ∇ 2 φ t+∆t i,j -(φ t y )i,j (φ t x )i,j
, where we used usual finite differences for the Hessian matrix ∇ 2 φ and a first order convex scheme [START_REF] Sethian | Level set methods[END_REF] for the advective term ∇φ • w. This scheme enables to compute the surface gradients in the appropriate direction: ∇φ T i,j wi,j = max(ui,j, 0)(φi,j) - x + min(ui,j, 0)(φi,j) + x + max(vi,j, 0)(φi,j) - y + min(vi,j, 0)(φi,j) + y ,

where (φ) - x and (φ) - y are the left semi-finite differences, whereas (φ) +

x and (φ) + y are the right ones. The discrete linear tangent operator is similarly defined as:

∂ φ t i,j Mδφ t+∆t i,j = M φ t i,j δφ t+∆t i,j - 2ε (A B) ||∇φ|| 4 (δφ t+∆t x )i,j (δφ t+∆t y )i,j ,
where A and B are: As for the iterative solver involved in the implicit discretization, we used a conjugated gradient optimization. The discretization of the adjoint evolution model is finally obtained as the transposed matrix corresponding to the discretization of the derivative of the evolution law operator.

A = φ t x φ t

Initial condition

In order to define an initial condition, we assume that an initial contour of the target object is available. Such contours may be provided by any segmentation technique or specified by the user. Given this initial contour, we assigned a signed distance function to the implicit function at the initial time. This initial condition is also assumed to hold up to some uncertainty.

More precisely, the value of φ(x, t 0 ) is set to the distance g(x, Γ(t 0 )) of the closest point of a given initial curve Γ(t 0 ), with the convention that g(x, Γ(t 0 )) is negative inside the contour, and positive outside. An additive control variable models the uncertainty we have on the initial curve. This initial condition reads: φ(x, t 0 ) = g(x, Γ(t 0 )) + η(x, t).

The covariance matrix B models the uncertainty on the initial state has been defined as a diagonal matrix B(x) = Ide -|g(x,t0)| . This covariance fixes a low uncertainty on the level set in the vicinity of the initial given curve. This uncertainty increases as soon as we move away from the initial contour.

Measurement equation

Associated to the evolution law and to the initialization process we previously described, we have to define a measurement equation linking the unknown state variable and available observations. In this work, the observation function Y (t) is assumed to be related to the unknown state function φ through a differential operator H and an additive discrepancy function ǫ:

Y (x, t) = H(φ(x, t)) + ǫ(x, t). (35)
In this work we use two different measurements equation for the curve tracking issue. Let us further describe the two options we devised.

Noisy curves observation

The simplest measurements equation that can be settled consists to define an observation function, Y (t), in the same space as the state variable φ, with the identity as measurements

Papadakis & Mémin

operator, H = Id. Sticking first to this case, we define the observation function as the signed distance map to an observed noisy closed curve, g(x, Γ(t)). These curves Γ(t) are assumed to be provided by a basic threshold segmentation method or by some moving object detection method. Let us outline that such observations are generally of bad quality. As a matter of fact, thresholding may lead to very noisy curves and motion detectors are usually prone to miss-detection or may even sometime fails to provide any detection mask when the object motion is too low. In this case, we define the covariance R of the measurements discrepancy as

R(x, t) = R min + (R max -R min )(Id -e -|Y (x,t)| ).
This function has lower values in the vicinity of the observed curves (R min = 10) and higher values faraway from them (R max = 50). When there is no observed curve, all the corresponding values of this covariance function are set to infinity.

Image observation based on local probability density

The previous observations require the use of an external segmentation process which we would like to avoid. To directly rely on the image data, a more complex measurement equation can

F (φ, I)(x, t) = [1 -d B (ρ Vx , ρ o )] 2 1 φ(x)<0 + [1 -d B (ρ Vx , ρ b )] 2 1 φ(x)≥0 = ǫ(x, t), (36)
where d B is the Bhattacharya probability density distance measure and ρ Vx is the probability density in the neighborhood of pixel x and ǫ is the function modelling the observation error.

By replacing the densities with intensity average, we retrieve the Chan and Vese functional proposed for image segmentation [START_REF] Chan | Active contours whitout edges[END_REF]. The Bhattacharya distance is in our case defined as:

d B (ρ 1 , ρ 2 ) = 255 0 ρ 1 (z)ρ 2 (z)dz.
This distance equals one if the two densities are identical. Let us note that the discrepancy between the measurements and the state variable does not appear explicitly. The measurement operator gathers here the measurements (the local intensity distribution) and the relation between the measurements and the state variable (the appropriate Bhattacharya distance).

The corresponding linear tangent operator in the sense of distributions [START_REF] Chan | Active contours whitout edges[END_REF] is:

∂ φ F = [1 -d B (ρ Vx , ρ o )] 2 -[1 -d B (ρ Vx , ρ b )] 2 δ(φ), (37)
where δ(•) is the Dirac function. In this case the covariance associated to the measurements error ǫ has been fixed to a diagonal matrix corresponding to the minimal empiric local

INRIA photometric covariance: R(x, t) = E[M in((1 -d B (ρ Vx , ρ o )) 2 , (1 -d B (ρ Vx , ρ b )) 2 )].
Finally, let us recall that for these first series of experiments the motion field transporting the curve of interest is assumed to be provided through an external motion estimation technique [START_REF] Mémin | Dense estimation and object-based segmentation of the optical flow with robust techniques[END_REF].

Experimental results

The first example shown on Figure 2 concerns a meteorological image sequence of the Meteosat infra red channel. We rely here on the simplest measurement model built from noisy curves. It is extracted through a simple threshold technique that select a photometric Infrared level line. Thus, we aim at tracking an iso-temperature curve. To demonstrate the robustness of our technique, we assume that this observation is missing at some instant of the sequence. The results and the observed curves are given in figure 2. The obtained results illustrate that the proposed technique keeps the adaptive topology property of level set methods, and in the same time, despite missing measurements, maintains a temporal consistency of the curves extracted. We then applied the technique to the tracking of a running tiger. This noisy sequence is composed of 27 frames and exhibits motion blur at some places. The measurements are provided here by local photometric histograms. The initial curve that determines probability density functions of the tiger and the background has been obtained with a simple threshold technique. The results obtained are shown on figure 5. For this sequence we also plotted on figure 6 a set of segmentation obtained through a Chan and Vese segmentation process [START_REF] Chan | Active contours whitout edges[END_REF] based on the same data model than our measurements model (eq. 36-37) together with an additional penalization term on the curve length. As can be observed from figure 6 the masks obtained with this spatial segmentation technique are of good quality for some frames.

t
t = 0 t = 2 t = 4 t = 7 t = 10 t = 12 t = 14 t = 17 t = 20 t = 22 t = 24 t = 26
They nevertheless appear to be unstable along time and would require a delicate tuning of the parameter to obtained a consistent sequence of curves. At the opposite, the curves provided by the proposed technique are more stable in time and consistent with respect to the object shape and its deformation. Compared to traditional segmentation techniques the assimilation techniques provides results which reflect in a more coherent way the topology and the deformation of the target object along time. Due to the bad quality of the image sequence and to the corresponding velocity fields, we can observe that it is nevertheless difficult for the curve to fit precisely and in a continuous way to the photometric boundaries of the object.

To further illustrate the differences between the results obtained through assimilation and successive photometric segmentations we finally run our method on a cardiac magnetic resonance imaging sequence the same photometric distribution (see figure 8). Incoherent merging or splitting of regions regarding the effective deformations of an object shape is maybe one of the main problems met when running spatio-temporal analysis on the basis of consecutive single spatial analysis (even chained together through their initializations). As for the computational cost of the method, it appears to be similar to standard level set segmentation. For instance, for a 512 image sequence of 20 frames the technique takes around two minutes of computation on a standard PC.

All the results shown here have been obtained considering as reliable external inputs a sequences of velocity fields describing the motion of a region of interest. In complex cases involving for instance occlusions, or situations for which the motion measurements are erroneous, this first assimilation system would not give satisfying results. During an occlusion, it is mandatory that the unobserved curve's motion fields are inferred on the basis of some evolution law. In other words, the curve's motion field must also be tracked. To that end, an extended state vector gathering both a motion field and a curve representation we be must considered. These coupled variables have to be assimilated. This is the topics of the following part of the paper.

Joint assimilation of motions and curves

The goal of this section is to show how the previous assimilation system can be enriched to handle the coupled tracking of a target curve and its motion. The idea consists in augment- ing the state variable with a motion representation. To that end, an appropriate motion conservation law and motion measurements must also be added to the overall system. Two very different cases can be distinguished: fluid motions and natural objects motions. As a

Initialization t = 1 t = 2 t = 3 t = 4 t = 5 t = 6 t = 7

Papadakis & Mémin

matter of fact, fluid motions are ruled by known and precise dynamical law and concerns situations for which the tracked curve does not delineate a specific object with its own motion.

Fluid motion is a global entity that is almost independent from the curve location. This case corresponds therefore to a weakly coupled assimilation problem. The case of motions of natural objects in video is completely different. In this case, trajectories may changed abruptly and the object's dynamics cannot be described by generic physical models. As the scene background and the object to track have their own motions, the curve delineating the target object is usually the place of strong motion discontinuities. A good estimation of the object motion or of the background motion is -at least in the vicinity of the object frontieressential. The quality of such estimation highly depends on the curve location. In this second case we thus have to tackle a tightly coupled assimilation problem.

These two cases will be studied separately. In the following section, we will first focus on the fluid motion case.

Fluid motion images

In several domains, the analysis of image sequences involving fluid phenomenon is of the highest importance. Such analysis is particularly meaningful as image sensors have the advantage to be non invasive and provide almost dense spatio-temporal observation of the flow. Satellite images of the atmosphere and oceans are of particular interest in geophysical sciences such as meteorology and oceanography where one wants to track cloud systems, estimate ocean streams or monitor the drift of passive entities such as icebergs or pollutant sheets. The analysis of fluid flow images is also crucial in several domain of experimental fluid mechanics for the analysis of peculiar experimental flows or for flow control purpose. For all the kinds of aforementioned applications and domains, it is of major interest to track features transported by the flow along time. Such a tracking which remains to estimate Lagrangian drifters for the structures of interest may be obtained from deterministic integration methods such as the Euler method or the Runge and Kutta integration technique. These numerical integration approaches rely on a continuous spatio-temporal vector field description and thus require the use of interpolation scheme over the whole spatial and temporal domain of interest.

As a consequence, they are quite sensitive to local errors measurements or to inaccurate motion estimates. When the images are noisy or if the flow velocities are of high magnitude and chaotic as in the case of turbulent flows, motion estimation tends to be quite difficult and prone to errors. Another source of error is inherent to motion estimation techniques (see for instance [START_REF] Barron | Performance of optical flow techniques[END_REF] for an extended review on motion estimation techniques). As a matter of fact, most of the motion estimation approaches use only a small set of images (usually two consecutive images of a sequence) and thus may suffer from a temporal inconsistency from frame to frame. The extension of spatial regularizers to spatio-temporal regularizers [START_REF] Bruhn | Lucas/kanade meets horn/schunck: Combining local and global optic flow methods[END_REF][START_REF] Papenberg | Highly accurate optic flow computation with theoretically justified warping[END_REF][START_REF] Weickert | Variational optic-flow computation with a spatio-temporal smoothness constraint[END_REF] or the introduction of simple dynamical constraint in motion segmentation techniques INRIA relies mainly on crude stationary assumptions. Dedicated fluid motion estimators introduce adequate data terms [START_REF] Corpetti | Dense estimation of fluid flows[END_REF], relevant basis functions [START_REF] Cuzol | A low dimensional fluid motion estimator[END_REF][START_REF] Cuzol | Vortex and source particles for fluid motion estimation[END_REF] or higher order regularizers [START_REF] Corpetti | Dense estimation of fluid flows[END_REF][START_REF] Kohlberger | Variational dense motion estimation using the helmholtz decomposition[END_REF][START_REF] Ruhnau | Variational optical flow estimation for particle image velocimetry[END_REF], in order to accurately estimate motion flows. However, these approaches do not include any dynamical constraint encouraging the preservation of conservation laws. Nothing warrants the solutions from being inconsistent with respect to the physical dynamical laws that the flow must obey.

As we will show it in the following section, the inclusion within the variational assimilation framework of some conservation laws combined with crude motion observation enables to improve notably the quality of the recovered dense motion fields.

Basic definitions and relations

ξ = ∇ ⊥ • w = ∂v ∂x - ∂u ∂y , ζ = ∇ • w = ∂u ∂x + ∂v ∂y . ( 38
)
The vorticity is related to the presence of a rotating motion, whereas the divergence is related to the presence of sinks and sources in a flow. Assuming that w vanishes at infinity2 , the vector field is decomposed using the orthogonal Helmholtz decomposition, as a sum of two potential functions gradient:

w = ∇ ⊥ Ψ + ∇Φ. ( 39
)
The stream function Ψ and the velocity potential Φ respectively correspond to the solenoidal and the irrotational part of the vector field w. They are linked to the divergence and vorticity maps through two Poisson Equations: ξ = ∆Ψ, ζ = ∆Φ. Expressing the solution of both equations as a convolution product with the Green kernel associated to the 2D Laplacian operator:

Ψ = G * ξ Φ = G * ζ, (40)
the whole velocity field can be recovered knowing its divergence and vorticity:

w = ∇ ⊥ G * ξ + ∇G * ζ △ = H G (ξ, ζ). ( 41
)
This computation can be very efficiently implemented in the Fourier domain.

Papadakis & Mémin

Fluid motion evolution equation In order to consider a tractable expression of the Navier-Stokes equation for the tracking problem, we rely in this work on the incompressible vorticity-velocity formulation of the Navier-Stokes equation:

∂ t ξ + w • ∇ξ = ν ξ ∆ξ. ( 42
)
This formulation states that the vorticity is transported by the velocity field and diffuses along time. For the divergence map, it is more difficult to exhibit any conservation law. We will assume that it behave like a noise. More precisely we will assume that the divergence map is a function of a Gaussian random variable, X t , with stationary increments (a Brownian motion) starting at points, x, of the grid where the initial velocity field admits a non null divergence. It can be shown through Ito formula and Kolmogorov's backward equation, that the expectation at time t of such a function, u(t, x) = E[ζ(X t )|X 0 ] obeys to a heat equation [START_REF]Stochastic differential equations[END_REF]:

∂ t u -ν ζ ∆u = 0, u(0, x) = ζ(x). (43)
According to this equation, we indeed make the assumption that the expectation of the divergence at any time of the sequence is a solution of a heat equation (i.e. it can be recovered from a smoothing of the initial motion field divergence map with a Gaussian function of standard deviation 2 √ ν ζ t ). The curl and divergence maps completely determine the underlying velocity field. Combining the motion equations ( 42) and ( 43) with the curve evolution law [START_REF] Martin | Influence of the noise model on level set active contours segmentation[END_REF] proposed in the previous section leads to a complete dynamical model for the state function X = [φ, ξ, ζ] gathering motions and curve descriptors, up to a control function ν:

∂ t     φ ξ ζ     +      ∇φ • w -ε ∆φ -∇ T φ∇ 2 φ∇φ ||∇φ|| 2 w • ∇ξ -ν ξ ∆ξ -ν ζ ∆ζ      = ν. (44)
The first component of this evolution law is discretized in the same way as before. The vorticity-velocity equation is discretized following a total variation diminishing (TVD) scheme proposed in [START_REF] Kurganov | A third-order semidiscrete central scheme for conservation laws and convection-diffusion equations[END_REF][START_REF] Kurganov | New high-resolution central schemes for nonlinear conservation laws and convection-diffusion equations[END_REF] and associated to a third order Runge-Kutta integration [START_REF] Shu | Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws[END_REF]. This discretization is described in appendix A. The divergence equation is integrated through a stable implicit scheme. The motion field is updated at each time step in the Fourier domain with equation [START_REF] Papenberg | Highly accurate optic flow computation with theoretically justified warping[END_REF]. Through this whole non-oscillatory scheme, the vorticity-velocity and divergence equations can be efficiently simulated.

Fluid flows observations

As fluid motions are assumed to be continuous and defined over the whole image domain, we will assume that an observation motion field w obs is available.

This motion field can be provided by any dense motion estimator. In this work, a dense χ, its linear tangent operator is itself. Nevertheless, as this operator is expressed through a convolution product, the expression of its adjoint ∂ χ H * G is not trivial. Following Appendix B it is defined as:

∂ χ H * G (w) = -H G (w).
The different covariance matrices corresponding to the motion components and involved in the dynamic, in the initialization and in the observation equations have been fixed to diagonal matrices with constant values (typically B = 0.1 and R = 10). As the dynamics is assumed to be quite accurate, the model covariance has been fixed to a low value of Q = 0.005 for all the fluid image sequence. The computational cost of this process is smaller than the optical flow computation on the whole sequence and it can be seen as a post-processing of fluid dense motions. Otherwise, as the system we settled for the fluid motion case is weakly coupled, the assimilation of the curve and motion components do not have to be done simultaneously. Since the motion components do not depend on the curve components, they are assimilated first. Once the sequence of velocity fields have been analyzed, the curve are in turn assimilated.

Natural object tracking

As explained in the introduction of the previous section, the situation of natural objects either rigid or deformable is completely different. Hence, the curve and its motion are tightly coupled. We have here to deal with an interlaced problem which can only be solved in a joint way.

No universal physical law can be stated for natural objects observed from videos. Therefore, it is always a great difficulty in tracking applications to design a general law describing their evolution with accuracy. General rough models such as constant velocity or constant acceleration dynamic models are commonly used [START_REF] Jackson | Tracking deformable moving objects under severe occlusions[END_REF]. This work also relies on that type of models. We consider two kinds of motion representation: a dense motion model related to the whole image and an affine velocity directly linked to the trajectory of the tracked object.

Dense motion

Evolution law for dense motion field We will assume that the targeted object move with a constant velocity along its trajectory, up to an additive control function. In differential form this reads:

dw dt = ν w . (45)
As w = [u, v] T , the evolution law for the dense motion fields is:

∂ t u + ∇u • w = ν u , ∂ t v + ∇v • w = ν v , (46)
where ν u and ν v are control function whose associated covariances are set to constant diagonal matrices. As the dynamic is here very likely approximative we used on the different sequence involved a value of Q u = Q v = 0.5 for these coefficients. Combining the previous evolution law of the curve component [START_REF] Martin | Influence of the noise model on level set active contours segmentation[END_REF] and the motion dynamic ( 46) leads to a coupled evolution model for the state function [φ, u, v] T :

∂ t     φ u v     +      ∇φ • w -ε ∆φ -∇ T φ∇ 2 φ∇φ ||∇φ|| 2 ∇u • w ∇v • w      = ν. (47)
To estimate the whole trajectory of the velocity field, it is necessary to define appropriate motion measurements.

Dense motion observations Accurate dense motion observations may be provided by optical flow technique at each frame of the sequence. Nevertheless, in case of occlusions, these motion measurements do not describe the tracked object trajectory and must not be taken into account. We describe in the following the mechanism on which we resort to properly assimilate a complete trajectory of velocity fields. Assuming that the whole object is visible on the two first frames, the initial velocity field estimated between these two frames constitutes a good representation of the initial object's motion. In the following frames, in order to only consider velocity measurements related to the object of interest, it is necessary to define informative values for the covariance matrix associated to the motion measurements.

High confidence values must be set at points where the object is visible (and low values where it is hidden). To fix these values, two distinct cases depending on the kind of contours observation we are dealing with must be considered separately.

-Noisy curve observations model: In such a measurements model (3.3.1), noisy or possibly inaccurate curves are assumed to be available when the whole object or some part of it is visible (i.e Y obs (x, t) < 0). This provides area of the image where the object is visible. The motion observation covariance matrices can then be defined through the level set surfaces corresponding to these observed curves:

R(x, t) = R min + (R max -R min )1 Y obs (x,t)≥0 .
Hence, the dense motion observations are preponderantly considered only in regions where the object is visible. In our applications, we systematically chose R min = 10 and R max = 50.

-Image observation model: When the measurements model relies on image local histograms (3.3.2), the situation is more problematic as no curves are directly observed.

We also want to avoid the use of an external segmentation process to infer occlusion situations. Considering the object and the background reference histograms (ρ o and ρ b ), it is possible to define adequate confidence values for the motion observation covariance matrix. We propose to fix them as:

R(x, t) = R min + (R max -R min )1 dB (ρV x ,ρo)<dB (ρV x ,ρ b ) . (3.3.2) (48)
Thus, low values (R min = 5) are associated to points that likely belong to the background, whereas higher confidence values (R max = 20) are related to the other points.

Affine velocity

In some cases, it is sometimes sufficient to rely on a reduced parametric representation of the object's motion. We now consider an affine velocity directly related to the target. In the following, we will assume that the apparent motion over the area delineated by the curve follows an affine parametric motion, i.e:

w(x, t) = S(t)x + T (t) ∀x∈ Ω|φ(x) < 0, [START_REF] Weickert | Variational optic-flow computation with a spatio-temporal smoothness constraint[END_REF] where S(t) is the rotation, divergence and shear matrix and T (t) is a translation vector.

Evolution law for affine velocity Assuming that the targeted object move with a constant velocity along its trajectory, up to a control function, as defined in equation ( 45), we obtain:

∂ t S(t) + S 2 (t) x + ∂ t T (t) + S(t)T (t) = ν w ∀x ∈ Ω|φ(x) < 0. ( 50
)
As this equation must be respected for all points inside the object, we get the following system of equations:

∂ t S(t) + S 2 (t) = ν S , ∂ t T (t) + S(t)T (t) = ν T , (51)
where ν S and ν T are control functions with constant diagonal covariance matrices (we typically chose Q S = Q T = 0.01). If S -1 (t 0 ) exists, then a closed form solution of the system of equations (51) exists:

S(t) = (S -1 (t 0 ) -Id t) -1 T (t) = S(t)T -1 (t 0 )T (t 0 ) (52)
Papadakis & Mémin

These solutions can be readily used for the forward integration of the motion components in the assimilation process. Combining the curve evolution law [START_REF] Martin | Influence of the noise model on level set active contours segmentation[END_REF] and the previous motion dynamical law (51) leads to a coupled evolution model for the state function [φ, S, T ] T :

∂t     φ S T     +      ∇φ • w -ε ∆φ -∇ T φ∇ 2 φ∇φ ||∇φ|| 2 S 2 (t) S(t)T (t)      = ν. ( 53
)
Let us explain now the different measurements Y S and Y T associated to this evolution law.

Affine velocity observations

As for the dense motion model, it is essential to have some significant motion measurements at one's disposal. Both measurement models should still be distinguished.

- For such a measurement model, we have specifically considered a final condition on the object contour as an additional equation of the assimilation system (see section 2.3). Assuming that the object is visible (at least partially) both on the initial and final image of the sequence, motion measurements associated to these two object locations allow us to assimilate roughly as a first guess the object motion. Latter on, the current sequence of curve estimates provides support to run a robust parametric estimation and thus enables to refine the motion measurements. More precisely, the motion estimation is run only on points of the support region associated to high confidence covariance value. The covariance values are fixed in the same way as defined in [START_REF] Talagrand | Variational assimilation of meteorological observations with the adjoint vorticity equation[END_REF]. Such a recipe avoids to perform motion estimation on occluded parts of the target.

Let us remark that the computational cost of the affine velocity assimilation is very small in comparison with the dense motion case. Indeed, the motion state variables are only (2 × 2) matrix and (1 × 2) vector, so the assimilation process for affine motion is instantaneous. As for the computational cost of a dense velocity field assimilation it remains smaller than the computation of the optical flow on the whole image sequence. It can be observed that the assimilation not only denoises the observations, but also allows to recover some small scales structures thanks to the vorticity-velocity dynamical model. The maps, these measurements constitute very noisy motion observations for the system we want to track. As we used a fast motion estimator, the motion observations present a lot of temporal inconsistencies with non plausible discontinuities artifacts. The photometric levellines are also quite noisy. To further demonstrate the ability of the proposed method, these photometric level lines have not been extracted for each images (these level lines have been extracted only every three images along the sequence). The estimated curve superimposed on the initial image sequence and the recovered vorticity maps are respectively shown on columns (c) and (d). We can notice that the motion fields recovered after assimilation are much more coherent from a physical point of view.

Results

The results presented in figure 13 shows the final curves recovered for the Vince cyclone.

This sequence is composed of 19 images obtained the 9 th of October 2005 from 00H15 to 4H45.We thank the Laboratoire de Météorologie Dynamique for providing us this sequence.

For this image sequence, we rely on the local photometric distributions measurement model.

That is to say we aim at tracking here a temperature distribution of reference defined from an initial mask on the first image. The joint assimilation of motion and curve enables to track perfectly the cyclone. Natural objects To demonstrate the potential benefit of a joint assimilation of closed curves and their motions for natural video sequences, we show results on two sequences exhibiting long term total occlusions. As a first result we show the tracking of two curves representing the photometric level lines of two cars in a video sequence of 27 frames. As can be observed, one of the car is occluded by the other during a long time period. In this case, we rely on an affine parameterization to describe the the objects motion. The motion measurements are given by motion estimates within the areas delineated by the observed curves when they are available. The level lines corresponding to the curve measurements are shown on the column (a) of the figure 14.

Obviously, there is neither curve nor motion measurements corresponding to the white car during its occlusion. Two level functions (one for each object) have been used in this example.

Column (b) of figure 14 present the curves recovered for the two cars. As can be observed these level lines are of bad quality and do not really described a precise boundary of the objects of interest. Let us note that even if the white car does not have a constant velocity during the occlusion, the assimilation process enables to predict quite well its photometric level line. For the tracking of natural objects in video, the main trouble concerns the velocities. If no coherent motion can be computed on the whole sequence as a first guess, the tracking algorithm will fail. Actually, a solution could be to rely on particle filtering, in order to find a first coherent motion. In fluid motion case, the proposed technique gives very good results.

It enables to enforce the motion field to respect some physical conservation and improve consequently the consistency of the solution.

Concerning applications, we believe it is of major interest when one aims at analyzing the deformations along time of entities enclosed by a curve. For instance, such an analysis is informative in environmental sciences for analyzing, representing or understanding complex phenomena. These phenomena may be only partially observable and subject to occlusions.

In video analysis, such a framework could be interesting for video post-processing of moving objects, such as object removal, object based video editing, off-line structuring, indexing of video databases.

Acknowledgments.

This work was supported by the European Community through the IST FET Open FLUID project (http://fluid.irisa.fr).

A Vorticity-velocity formulation discretization

Discretizing the vorticity-velocity equation ( 42) requires a careful handling. In particular, the advective term ∇ξ • w must be treated cautiously. A lot of non-oscillatory schemes based on conservation laws have been developed to solve this problem [START_REF] Harten | Uniformly high order accurate essentially non-oscillatory schemes, 111[END_REF][START_REF] Kurganov | New high-resolution central schemes for nonlinear conservation laws and convection-diffusion equations[END_REF][START_REF] Nessyahu | Non-oscillatory central differencing for hyperbolic conservation laws[END_REF]. Such schemes aim at considering the cell average and discretizing the value of the advective term at the boundaries of the cell with more accuracy than do classical finite difference schemes. Indeed, the gradients of the transported quantity are computed with high orders reconstructions in order to limit the slops and prevent from inappropriate numerical amplifications. The ENO (Essentially non-oscillatory) or WENO (Weighted ENO) are the most used schemes for reconstruction [START_REF] Levy | A third order central weno scheme for 2D conservation laws[END_REF][START_REF] Shu | Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws[END_REF]. To discretize the advective scheme, the velocity must be considered according to its local directions. The conservation laws are respected with the use of monotone flux; wellknown methods are Godunov and Lax-Friedricks flux. These techniques enable the scheme to be Total Variation Diminishing (TVD). It means that the total energy of the transported quantity can not increase along time. Hence, the TVD property prevents from numerical explosions and also assures that the shocks will be preserved. All these methods are well detailed in [START_REF] Shu | Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws[END_REF]. In our work, the reconstruction of the vorticity on the cell boundaries is realized through a second order method [START_REF] Levy | Non-oscillatory central schemes for the incompressible 2-d euler equations[END_REF] based on a Min-Mod limiter: are computed in the same way. As the Mid-Mod limiter choose the smallest slope, the reconstructed values of the vorticity on the cell boundaries limit the variation of the spatial gradients. To deal with the advective term, we use the following semidiscrete central scheme [START_REF] Kurganov | A third-order semidiscrete central scheme for conservation laws and convection-diffusion equations[END_REF][START_REF] Kurganov | New high-resolution central schemes for nonlinear conservation laws and convection-diffusion equations[END_REF]

with a numerical convection flux derived from the monotone Lax-Friedricks flux: This second order, semidiscrete central scheme is weakly dissipative and follows a TVD property [START_REF] Kurganov | A third-order semidiscrete central scheme for conservation laws and convection-diffusion equations[END_REF][START_REF] Levy | Non-oscillatory central schemes for the incompressible 2-d euler equations[END_REF]. The intermediate values of the velocities are computed with a fourth-order averaging:

H x i+
u i+ 1 2 ,j (t) = -ui+2,j(t) + 9ui+1,j (t) + 9ui,j (t)ui-1,j (t) 16

v i,j+ 1 2
(t) = -vi,j+2(t) + 9vi,j+1(t) + 9vi,j(t)vi,j-1(t) 16 .

(56)

The linear viscosity ∆ξ is approximated by the fourth-order central differencing:

Di,j(t) = -ξi+2,j(t) + 16ξi+1,j (t) -30ξi,j(t) + 16ξi-1,j (t)ξi-2,j (t) 12∆x 2 + -ξi,j+2(t) + 16ξi,j+1(t) -30ξi,j(t) + 16ξi,j-1(t)ξi,j-2(t) 12∆y 2 .

(57)

The time integration is realized with a third-order Runge Kutta scheme, which also respect the TVD property [START_REF] Shu | Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws[END_REF].

B Fluid motion adjoint observation operator

In order to have more compact notations, let Z(x) = [Ψ(x), Φ(x)] T denotes the potential coordinates. These coordinates are linked to the div-curl coordinates, χ(x) = [ξ(x), ζ(x)] T through relation [START_REF] Papadakis | Variational optimal control technique for the tracking of deformable objects[END_REF]:

Z = G * χ. (58)
Besides, relation [START_REF] Papenberg | Highly accurate optic flow computation with theoretically justified warping[END_REF] reads immediately:

u = ∇ ⊥ • Z and v = ∇ • Z. (59)
Similarly to what has been done before, the vector field Z may be written on the basis of the Helmholtz decomposition as the sum of the gradient of two scalar functions. Relation (59) allows us to write both motion components as solutions of two Poisson equations. Expressing these solutions through Green kernel enables to write the vector field Z as:

Z = ∇ ⊥ G * u + ∇G * v = H G (u, v). (60)
It can be can noticed that the same operator is involved in [START_REF] Papenberg | Highly accurate optic flow computation with theoretically justified warping[END_REF] and (60). These last notations, which link the stream function and the velocity potential to the vector field, will be useful in the following. Proof:

w A , w B (59) = u A v A • ∇ ⊥ • Z B ∇ • Z B = u A (Φ B x -Ψ B y ) + v A (Φ B y + Ψ B x ) (Integration by parts) = -Ψ B (v A x -u A y ) -Φ B (u A x + v A y ) + Boundary terms = 0 = -χ A , Z B
Using this property, the adjoint of the linear tangent operator can be found:

w A , w B (41) = H G (χ A ), w B = ∂ χ H G (χ A ), w B def = χ A , ∂ χ H * G (w B ) w A , w B P r = -χ A , Z B

  Méthode variationnelle pour le suivi couplé de courbes et de mouvement Résumé : Dans cet article, nous décrivons un nouveau procédé permettant de suivre des courbes fermées et leur champ de mouvement associé. La méthode proposée réalise le suivi continu d'une courbe déformable et d'un champ de vitesse au cours d'une séquence d'images. Le problème consiste à minimiser une fonctionnelle spatio-temporelle globale par rapport à un ensemble de variables représentant la courbe et son mouvement associé. En s'appuyant sur le principe du contrôle optimal, une loi d'évolution est intégrée le long de la séquence, puis sa loi d'évolution adjointe est intégrée rétrogradement. Cette dernière équation aux dérivées partielles est liée à un terme de distance entre les variables d'état fournies par la première intégration et des mesures discrètes du système. Les courbes fermées sont représentées par des surfaces implicites. Le mouvement est quant à lui décrit, suivant l'application, par un champ de vecteurs ou par des cartes de vorticité et de divergence. L'efficacité de l'approche est démontrée sur deux types de séquences d'images présentant des objets déformables et des mouvements fluides.

Figure 1 :

 1 Figure 1: Assimilation algorithm principle. This figure gives a synoptic of the overall principle of the method. After an integration of the initial condition X0 (blue trajectory), a backward integration

  be built to deal with local probability distributions of the intensity function. This model compares at each point of the image domain Ω a local photometric histogram to two global probability density functions ρ o and ρ b modelling respectively the object and background intensity distribution. These two distributions are assumed to be estimated from the initial location of the target object. The measurements equation we propose in this case reads:

Figure 2 :Figure 3 :Figure 4 :

 234 Figure 2: Clouds sequence. Top row: Sample of observed curves. Bottom row: Recovered results superimposed on the corresponding image To demonstrate the interest of the second observation model based on the local probability density, we applied our assimilation technique with such measurements for the tracking of

Figure 5 :

 5 Figure 5: Tiger sequence. Result of the assimilation technique with the photometric measurement model based on local probability density.

1 .Figure 6 :

 16 Figure 6: Tiger sequence. Successive segmentations obtained through a Chan and Vese level-set techniques with a data model based on local probability density measurement and a Bhattacharya distance (eq. 36-37).

Figure 7 :Figure 8 :

 78 Figure 7: Cardiac magnetic resonance imaging sequence. Result of the assimilation technique with the photometric measurement model based on local probability density. The initial target region is shown in the first image of the top row.

  For fluid flows, the Navier-Stokes equation provides a universal general law for predicting the evolution of the flow. In this work, the formulation of the Navier-Stokes on which we will rely on, uses the vorticity ξ and the divergence ζ of a motion field w = [u, v] T :

INRIA

  motion field estimator dedicated to fluid flows is used[START_REF] Corpetti | Dense estimation of fluid flows[END_REF]. Referring to equation[START_REF] Papenberg | Highly accurate optic flow computation with theoretically justified warping[END_REF], the observation operator for the motion components H(w) is defined as w obs = H G (ζ, ξ). This operator links the motion measurements w obs to the state variable χ = [ζ, ξ]T (where χ is the compact notation for the div-curl coordinates). The observation equation for the motion components finally reads Y (χ) = H G (χ). The observation operator is linear w.r.t

-

  Noisy curve observations model: For the measurements model (3.3.1), the observed contours enable to perform a robust parametric motion estimation inside the region of interest. Such a process directly provides the components of the observation function, Y S and Y T . The values of the least square residuals are used to fix the corresponding covariance values [37]. When no curve are provided or if the area delineated is to small, no motion estimations are considered and the corresponding values of the observation covariance matrix are set to infinity. Image observation model: When the measurements model relies directly on image local histograms (3.3.2), the only region available on which it is possible to get motion estimates is the one specified at the initial time as the initial condition. This single rough parametric motion measurement reveals to be generally unable to anticipate sufficiently well the object's trajectory on the whole sequence during the first forward step.

Fluid motionFigure 9 :

 9 Figure 9: 2D Direct Numerical Simulation .a) Particle images sequence. b) True vorticity. c) Vorticity observed by optic flow estimator. d) Assimilated vorticity.

  been provided by a dedicated optic flow estimator[START_REF] Corpetti | Dense estimation of fluid flows[END_REF]. To give some quantitative evaluation results, we present in figure10comparative errors between the ground truth and the obtained results. On this example, the assimilation process enables to significantly improve the quality of the recovered motion field (about 30%).

Figure 10 :Figure 11 :Figure 12 :

 101112 Figure 10: Comparison of errors. The red curve outlines the mean square error of the vorticity computed by the optic flow technique on the 50 images of the turbulent particles image sequence.The green curve exhibits the error obtained at the end of the assimilation process. The assimilated vorticity is then closer to the reality than the observed one.

Figure 13 :

 13 photometric measurement model based on local probability density.

Figure 14 : 5 ConclusionFigure 15 :

 14515 Figure 14: Car sequence 2. (a) Sample of observed photometric level lines. (b) Curves after joint assimilation.

  infi(xi) if xi ≥ 0 ∀i sup i (xi) if xi ≤ 0 ∀i 0 otherwiseThe terms ξ + i,j+ 1 2 and ξ - i,j+ 1 2

Property B. 1

 1 Let w A and w B be two vector fields belonging to C 2 (Ω) ∩ L 2 (Ω), where Ω is a bounded set of R 2 with border ∂Ω. Assuming that the components of both vector fields admit Dirichlet boundaries conditions, i.e u A,B | ∂Ω = v A,B | ∂Ω = 0,then w A , w B = -Z A , χ B =χ A , Z B .

  by identifying the terms∂ χ H * G (w B ) = -Z B = -H G (w B). It finally reads:∂ χ H * G (w) = -H G (w).(61)

  :

	∂ξi,j ∂t	= -	H x i+ 1 2 ,j (t) -H x i-1 2 ,j (t) ∆x	-	H y i,j+ 1 2	(t) -H y i,j-1 2 ∆y	(t)	+ ν ξ Di,j

RR n 6283

Papadakis & Mémin

http://mrel.usc.edu/class/preview.html RR n 6283

In order to restrict the computation to the image domain, a divergence and curl free global transportation component is removed from the vector field. This field is estimated on the basis of a Horn and Schunck estimator associated to a high smoothing penalty[START_REF] Corpetti | Dense estimation of fluid flows[END_REF].RR n 6283