Christophe Avenel
email: christophe.avenel@lip6.fr

Etienne Mémin
email: etienne.memin@inria.fr

Patrick Pérez
email: patrick.perez@technicolor.com

Stochastic level set

Introduction

The inference of the trajectory along time of an interface delineating two regions is of major interest in various domains such as computer vision, medical imaging or fluid flow analysis. For the latter domain, when complex flows such as oceanic or atmospheric flows are involved, the forecasting of such an interface on the sole basis of a deterministic numerical model is extremely difficult due to the underlying physics. As a matter of fact, the corresponding physical models are in general only partially known and involve complex small scales phenomena that are difficult to model at large scales. In addition, when one aims to track within the image plane a 3D curve delineating two photometric areas of interest, the formalization of such an evolution model is an intricate issue. Even if it is conceptually possible, this task would indeed involve 3D numerical models at a scale of the same order as the pixel resolution. For most geophysical flows observed through satellite sensors with ever-increasing resolution capabilities, handling numerical models at such scales (i.e., in the range of the submeso scales that lie in between 1km and 50km) represents a formidable computational task. The establishment of empirical dynamical models learned from image data through appropriate data assimilation procedure constitutes hence an interesting alternative. It might be useful in numerous applications such as pollutant sheets tracking and forecasting, potentially dangerous cloud system monitoring or complex physical phenomena evolution analysis over long time periods for climatic variation studies. Such a tracking -or data assimilation -objective faces two major difficulties. First of all, as we aim at inferring from noisy image data the whole trajectory of a closed curve, it is necessary to define a technique that allows the robust tracking of such an interface along time. This tracking procedure defined from data assimilation strategy incorporates as a crucial ingredient an evolution model of the variable of interest. As we do not wish to rely on exact physical evolution models, but rather to consider a simplified dynamics with uncertainty terms, it is quite natural to consider a stochastic assimilation framework and hence a stochastic dynamics. The constitution of a sound stochastic dynamics for the interface constitutes the second difficulty of such a study.

As we do not want to focus in this study on the tracking of particular phenomena but instead to deal with regions potentially undergoing any kind of complex deformations including those accompagnied by topological changes, we will confine ourselves to a level set representation of the curve. This representation has the great advantage to formulate the curve evolution within an Eulerian framework. In particular, it handles easily topological changes as opposed to Lagrangian curve representations (e.g., spline curves) that require explicit split/merge mechanisms. Since the seminal work of Osher and Sethian [START_REF] Osher | Fronts propagating with curvature dependent speed: Algorithms based on Hamilton-Jacobi formulation[END_REF] , this representation, originally proposed by Dervieux and Thomasset [START_REF] Dervieux | A finite element method for the simulation of rayleightaylor instability[END_REF] for fluid interface simulation, has been an intensive subject of research in domains like computer graphics, computational fluid dynamics, or computer vision (see for instance the textbooks [START_REF] Osher | Level set methods and dynamic implicit surface[END_REF][START_REF] Osher | Geometric level set methods in Imaging, Vision and Graphics[END_REF][START_REF] Sethian | Level set methods[END_REF] and the references therein). This representation has motivated numerous developments in computational fluid dynamics and shape optimization to describe the motion of an interface [START_REF] Chang | A level set formulation of eulerian interface capturing methods for incompressible fluid flows[END_REF], the formalization of fluid structure interactions [START_REF] Coqrelle | A vortex level set method for the two-way interacting coupling of an incompressible fluid with colliding rigid bodies[END_REF][START_REF] Cottet | A level set method for fluid structure interactions with immersed surfaces[END_REF], or the study of elastic structures [START_REF] Allaire | Structural optimization using sensitivity analysis and a level-set method[END_REF]. In computer vision, level sets have been tremendously used for image segmentation [START_REF] Caselles | Geodesic active contours[END_REF][START_REF] Chan | Active contours whithout edges[END_REF][START_REF] Vese | A multiphase level set framework for image segmentation using the Mumford and Shah model[END_REF] while, surprisingly enough, they have been used only occasionally for tracking problems [START_REF] Kimmel | Tracking level sets by level sets: a method for solving the shape from shading problem[END_REF]. In addition, most of the approaches proposed so far for the tracking of a level set representation are simply defined as successive, almost independent, detection processes on each image of the sequence [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] Niethammer | Dynamic geodesic snakes for visual tracking[END_REF][START_REF] Paragios | Geodesic active regions: a new framework to deal with frame partition problems in computer vision[END_REF]. Even if those techniques might include eventually a temporal smoothing, they can't truly be considered as tracking processes as they do not guaranty any temporal coherence along the curve points trajectory. Such temporal incoherences are all the more pregnant since ambiguities due to clutter noise, missing data (occlusions) or illumination variations are observed. In addition to this, the few probabilistic techniques that have been proposed so far in the literature for curve tracking [START_REF] Dambreville | Tracking deformable objects with unscented kalman filtering and geometric active contours[END_REF][START_REF] Rathi | Tracking deforming objects using particle filtering for geometric active contours[END_REF] are built upon ad-hoc linear models of the curve evolution law, suited only to the tracking of objects undergoing small or quasi-rigid deformations. Such approximations are totally inappropriate for the highly nonlinear evolution laws driving the phenomena we focus on in this study. Let us note that, to alleviate most of the limitations mentioned above, a deterministic optimal control approach has been proposed in [START_REF] Papadakis | A variational technique for time consistent tracking of curves and motion[END_REF]. This level sets variational tracking approach constitutes nevertheless a batch approach of the curve tracking problem, while we propose here an online recursive technique enabling the analysis of long sequences.

In the approach we propose, the level set evolution law takes thus the form of a stochastic partial differential equation. A key ingredient is the noise term that introduces an anisotropic smoothing with a mean curvature motion component complemented with a Hessian diffusion. The drift term responsible of the transport of the curve's implicit surface will be learned from the past trajectories and the image data through a data assimilation procedure. The drift is expressed through an additional vectorial level sets maintaining along time the correspondences between the curve's points [START_REF] Pons | Maintaining the point correspondence in the level set framework[END_REF]. The dynamics describing the evolution of this auxiliary function incorporates the uncertainties on the curve deformations as well. The curve dynamics is supplemented with a local photometric information that takes the form of a data-driven force, in order to guide more efficiently the curve prediction toward the next image observation.

As for the data assimilation process, it is formulated as a stochastic filtering problem in which the available photometric data are filtered by the level set stochastic evolution law. For the reasons mentioned previously stochastic filtering constitutes a very appealing framework for data assimilation. In high dimensional spaces those filters are unfortunately notably difficult to implement efficiently with particle or ensemble filters. To mitigate this difficulty, it is very important to devise dynamics that are the most accurate as possible and in the same time to circumscribe the space of the random deformations applied on the curve. This indeed enables to efficiently draw samples but also makes possible an efficient exploration of the considered curve's state space. To that end, the curve dynamics on which we rely is formulated as a stochastic transport equation defined directly on the implicit surface function and includes constant displacement uncertainties along the curve tangent and normal directions. The method we propose is related to the technique proposed by Avenel et al. [START_REF] Avenel | Tracking closed curves with non-linear stochastic filters[END_REF]. However opposite to this work, the dynamics we propose does not rely on external velocity measurements. Furthermore our modelling is more general as we will show, since we do not implicitly introduce a mean curvature component to the transport velocity field.

The remainder of the paper is summarized as follows. We first introduce our new dynamical model, and the associated deterministic transportation drift. We then show how we implemented this dynamical model within a Bayesian particle filter, and what kind of measurements have been considered for this tracking problem. The last part of the paper shows various experimental results obtained with this technique. This paper corresponds to an extension of a conference papers [START_REF] Avenel | Tracking levels representation driven by a stochastic dynamics[END_REF], where only the general results were announced but not detailed. It also provides a more general derivation of earlier works [START_REF] Avenel | Tracking closed curves with non-linear stochastic filters[END_REF][START_REF] Avenel | Stochastic filtering of level sets for curve tracking[END_REF], where the level set surface was assumed to be time differentiable, which is not in general. This notion is of great importance since it impacts the differentiation formulas that must be used. We elaborate also on the relation between the proposed stochastic dynamics and the standard levels set evolution equation and provides results on very long sequences of daily satellite data with strong noise and even a repeated absence of data during summertime (hence over a quarter of the total duration).

Level set dynamics for closed curve evolution model

This section will first present the dynamical model we introduce, and then the method used to retrieve the velocity of the curve. Let us note that even if only the 2D curve's case is considered an extension to 3D surfaces can be however immediately derived from the very same principle.

Stochastic dynamics

In order to alleviate random sampling in high dimensional state spaces problem, we will restrain the potential uncertainty on the curve's deformations to belong to a space of low dimension. To that end, the uncertainty on the curve deformation will be defined through two independent constant Brownian motions directed along the curve's normal and tangent. The rationale is to keep all the dimensionality of the state variable but to add a perturbations of very small dimension in order to be able to sample it efficiently.

More precisely, we set for the curve, C t , a deformation model defined as:

dC t = w n ndt + σ n ndB n,t + σ τ n ⊥ dB τ,t . (1)
with dB n,t and dB τ,t two uncorrelated standard Brownian motions, σ n and σ τ two diffusion coefficients, n the unit vector normal to the curve, n ⊥ the unit vector tangent to the curve and the drift strength w n = w T n the projection on the curve's normal of a deterministic transport velocity field w. In this model, the random Brownian diffusion terms encode the uncertainties we have on the curve deformation. The first one is relative to a scale uncertainty whereas the second one, the curve tangential component, affects only the curve parametrization and is generally omitted in deterministic curve evolution model. In our case this term will have an influence on the final form of the corresponding level set evolution and hence must be kept. In addition, as we will define the drift transportation component, w n , from points curve correspondences, this term will be of importance and cannot be discarded. This drift component will be further detailed in the second part of this section. Let us remark, that non trivial tangential velocity may be also added into the geometric law (1) in order to improve and stabilize the numerical schemes of the interface evolution (see [START_REF] Hou | Removing the stiffness of curvature in computing 3D filaments[END_REF][START_REF] Hou | Removing the stiffness from interfacial flows and surface tension[END_REF][START_REF] Mikula | A direct method for solving an anisotropic mean curvature flow of plane curves with an external force[END_REF]). Let us note also, that the diffusion coefficients are not necessary constant along time. They may also eventually depends on the curve, nevertheless for sake of readability we will not indicate those dependencies in the notation unless specifically required.

In order to cope with curve's deformation of any kinds, it is essential to rely on a curve representation allowing handling topological changes that occur when the region of interest splits apart in several separated components or on the contrary reassembles. The level set formalism [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] has been specifically introduced in that goal to bypass the deficiency of splines based curves representation in such situations. In this representation, the curve C t at time t is defined as the zero level set of a scalar function ϕ(x, t) : Ω × R + → R:

C t = {x ∈ Ω | ϕ(x, t) = 0} , (2)
where Ω stands for the image spatial domain. The implicit surface function, ϕ, is chosen so as to have for instance positive values inside the curve and negative values outside. A common choice for the implicit function initialization is the signed distance function associated to initial curve but any other surface whose level set fits the initial curve of interest is possible. The surface evolution is then defined in order to stick at all times to the contour dynamics. This representation has the great advantage to allow describing, through a single implicit surface, a set of non-intersecting closed curves. The main geometric features of the curve that will be needed for the curve evolution computation can be directly obtained from the implicit surface. In particular, the inward unit normal and the mean curvature are respectively given by:

n = ∇ϕ |∇ϕ| and κ = ∇ • ∇ϕ |∇ϕ| = 1 |∇ϕ| (∆ϕ - 1 |∇ϕ| 2 ∇ϕ T ∇ 2 ϕ∇ϕ), (3)
where ∆ϕ and ∇ 2 ϕ denote the Laplacian and the Hessian of ϕ, respectively. Let us consider a C 2 (Ω) function ϕ : R 2 × R + → R, defined as a level set representation of the curve, which for any curve's point X t at time t is such that:

ϕ(X t , t) = 0. (4)
The surface ϕ can be used to express the deformation, dX t , of the curve (1) on the whole domain Ω:

dX t = w * n ∇ϕ |∇ϕ| dt + σ n ∇ϕ |∇ϕ| dB n,t + σ τ ∇ϕ ⊥ |∇ϕ| dB τ,t , (5)
where w * n is an extension to the whole image domain of the curve's drift strength and ∇ϕ ⊥ = (∂ y ϕ, -∂ x ϕ) T denotes the perpendicular gradient.

For a fixed point y ∈ Ω, the function ϕ is necessarily a random function of the form:

dϕ t (y) = b(y, t)dt + f (y, t)dB n,t + g(y, t)dB τ,t , (6)
where dϕ t (y) denotes the differential of ϕ(y, •) at a fixed grid point y. As a matter of fact, we can see that ϕ t (y) cannot be a deterministic function, since otherwise a standard Ito differentiation with a time scale separation between the rapid oscillating random terms and the slow deterministic components would yield a cancellation of the random term (σ n dB n,t = 0) and hence a null normal uncertainty, which is not the general goal followed here. Previous works based on Ito differentiation [START_REF] Avenel | Tracking closed curves with non-linear stochastic filters[END_REF][START_REF] Avenel | Stochastic filtering of level sets for curve tracking[END_REF] assumed thus implicitly only a noise directed along the curve tangent (σ n = 0) or a level sets conservation equation defined only up to ∇ϕσ n dB n,t , which implies a tracking of the zero level set only in mean. As ϕ(x, t) is a semi-martingale, the differential of ϕ(X t , t) has to be calculated through the Ito-Wentzell formula (see [START_REF] Kunita | Stochastic flows and stochastic differential equations[END_REF] theorem 3.3.1 on the differential of two stochastic processes composition ϕ • X t ). At point X t = x, we have:

dϕ(X t , t) = dϕ t (x) + ∇ϕ T dX + 1 2 i,j d X i t , X j t ∂ 2 ϕ ∂x i ∂x j dt + i d ∂ϕ ∂x i , X i t (7) = 0.
The term X i t , X j t denotes the joint quadratic variations of X i t and X j t , with X t = X 1 t , X 2 t . It is defined as the limit in probability over a partition {t 1 , . . . , t n } of [0, t] with 0 < t 1 < t 2 < • • • < t n < t, and a partition spacing

δt i = t i -t i-1 , noted as |δ n | = max i δt i and such that |δ n | → 0 when n → ∞: X i t , X j t = P lim |δ n |→0 t k ≤t (X i t k-1 -X i t k )(X j t k-1 -X j t k ), for n → ∞.
For Brownian motion, these quadratic covariations can be easily computed and are given by the following identities:

B i , B j = δ ij t a.s. h(t), h(t) = h(t), dB i = B j , h(t) = 0 a.s., (8)
where δ ij is the Kronecker symbol (δ ij = 1 if i = j, δ ij = 0 otherwise), and h(t) is a deterministic function.

At point X t = x we have thus ( 7) which reads:

dϕ(X t , t) = dϕ t (x) + |∇ϕ|w * n dt + σ n |∇ϕ|dB n,t + σ 2 n 2 1 |∇ϕ| 2 ∇ϕ T ∇ 2 ϕ∇ϕ dt + σ 2 τ 2 ∆ϕ - 1 |∇ϕ| 2 ∇ϕ T ∇ 2 ϕ∇ϕ dt + σ τ |∇ϕ| i ∂g ∂x i ∂ϕ ∂x i ⊥ dt + σ n |∇ϕ| i ∂f ∂x i ∂ϕ ∂x i dt (9)
= 0.

Let us note that compared to Ito differentiation the two last terms correspond to the cross variation of ∇ϕ and X t . These terms cancel for a deterministic function and we retrieve the traditional Ito formula. Let us remark that for correlated noise additional cross factors involving the product of the two diffusion coefficients would have to be included. Incorporating the expression (6) in equation ( 9), and equating respectively the random terms and the deterministic terms, enables immediately to get f (x, t) = -σ n |∇ϕ| and g(x, t) = 0. Now b(x, t) can be computed using

i ∂f ∂x i ∂ϕ ∂x i = -σ n |∇ϕ| ∇ϕ T ∇ 2 ϕ∇ϕ : b(x, t) = -|∇ϕ|w * n - σ 2 τ 2 ∆ϕ - 1 |∇ϕ| 2 ∇ϕ T ∇ 2 ϕ∇ϕ + σ 2 n 2 1 |∇ϕ| 2 ∇ϕ T ∇ϕ∇ϕ . (10)
And hence

dϕ t (x) = b(x, t)dt + f (x, t)dB n,t reads dϕ t (x) = -|∇ϕ|w * n dt - σ 2 τ |∇ϕ| 2 κdt + σ 2 n 2 1 |∇ϕ| 2 ∇ϕ T ∇ 2 ϕ∇ϕ dt -σ n |∇ϕ|dB n,t . (11)
We can remark that this stochastic partial differential equation involves a multiplicative noise, as the diffusion coefficient weighting the white noise term depends on ϕ. This level set evolution law can be simplified further if the motion field w * n is assumed to have a mean curvature flow component so as:

w * n = u * n -σ 2 τ |∇ϕ|κ. (12)
As it will subsequently explained in the next section, such an expression emerges immediately if the drift strength is defined as the gradient of a standard level set image segmentation functional. For such a drift strength we have then:

dϕ t (x) + |∇ϕ|u * n dt + σ n |∇ϕ|dB n,t = σ 2 τ dt 2 ∆ϕ + σ 2 n -σ 2 τ dt 2|∇ϕ| 2 ∇ϕ T ∇ 2 ϕ∇ϕ , (13)
which corresponds to the model considered in [START_REF] Avenel | Tracking closed curves with non-linear stochastic filters[END_REF]. For a null normal uncertainty, σ n = 0, we get a classical deterministic level set equation with a mean flow curvature motion component:

∂ t ϕ(x, t) + |∇ϕ|u * n = σ 2 τ 2 κ|∇ϕ|. (14)
For balanced uncertainties (σ n = σ τ ), the model corresponds to an advection diffusion equation with a stochastic (multiplicative) forcing:

dϕ t (x) + |∇ϕ|u * n dt = σ n dt 2 ∆ϕ -σ n |∇ϕ|dB n,t . (15)
Let us note that in all those cases the mean curvature flow component of the motion measurements does not have to be explicitly implemented. It could be imposed as an implicit assumption on the transportation field that is directly taken into account in the level set evolution. This enables a simple justification of the use of stochastic models that resemble very much in their formulation to classical level set evolution models. However, in the present work we rely on our more general stochastic representation (11) and we do not assume such an hypothesis for the transportation velocity fields. Nevertheless, as we will show later on, this mean curvature component could be very easily introduced through a data driven potential derived from a segmentation problem. The stochastic representation of the curve temporal evolution we derived in this section constitutes an essential ingredient of the tracking method proposed here as it will enable us to draw samples of forecasted deformed curves.

Drift component

The dynamics introduced in the previous section includes on the one hand, uncertainty terms that depends on the noise parameters (namely the diffusions coefficients) and on the other hand, an advection term that depends on a deterministic transportation field and eventually on the mean curvature as described in the previous section. In this section, we propose an expression of this transportation component that depends on the past curve's trajectory and also on the current image data in order to guide the random curve's samples toward meaningful areas of the state space.

The straightforward immediate idea consists in fixing the drift component from motion estimation techniques. This solution is often employed, but reveals to be less usefull than it might appear. It is for instance not at all adapted to handle missing data areas, where no motion measurements are available. Hence, areas of missing data, caused by occlusions or due to sensors failures in particular situations, require to set up detection mechanisms, which are not always that easy to define, but also to specify the motion missing values. In order to alleviate these difficulties, we propose instead to infer directly the curve's velocity field from each particle displacements, via a second implicit representation that keeps track of each curve point's initial location in the image plane. As proposed in [START_REF] Pons | Maintaining the point correspondence in the level set framework[END_REF], in order to implement this point correspondences, we introduce an additional vectorial level set ψ encoding on the spatial domain the curve's point location transportation in between two image frames. Keeping track of these backward point correspondences between the current evolving curve and a recent predecessor will allow us to derive an estimate of the curve's point velocity field. To that end, let us encode the Cartesian coordinates of the curve's points location at time k -1, through a vector-valued level set function ψ :

R 2 × R + → R 2 : ψ(x, k -1) = x. (16)
For time values within the interval between two consecutive images, [k -1, k], the vectorial function ψ(x, t) defines the location that point x ∈ Ω was occupying at the previous instant k -1. Function, ψ(x, t), is intrinsically attached to the curve and undergoes deformations imposed by the stochastic curve's evolution law [START_REF] Cottet | A level set method for fluid structure interactions with immersed surfaces[END_REF]. Its differential describing the transportation of the curve's coordinates at the previous instant yields at point X t = x a vectorial evolution equation of the same form as for the curve implicit surface:

dψ i (X t , t) = dψ i t (x) + (∇ψ i t ) T dX + 1 2 i,j d x i , x j ∂ 2 ψ i t ∂x i ∂x j + i d ∂ψ i t ∂x i , X i = 0.
(17) In the very same way as previously, identifying the terms of the semi-martingale associated to the differential dψ i t (x) at a fixed grid point through the consistency constraint [START_REF] Goldenberg | Fast geodesic active contours[END_REF], leads to:

dψ i t (x) = -(∇ψ i t ) T ∇ϕ |∇ϕ| w * n dt -(∇ψ i t ) T ( ∇ϕ |∇ϕ| σ n dB n,t + ∇ϕ ⊥ |∇ϕ| σ τ dB τ,t ) - A i dt 2|∇ϕ| 2 + σ n F i dt |∇ϕ| + σ τ G i dt |∇ϕ| . ( 18
)
The exact expression of terms A, F and G are detailed in the appendix. This vectorial evolution law enables us to define the curve's transportation component extended to the whole image plane, w * , (the asterisk notation indicates this extension of the curve drift) for the next inter frames interval [k, k + 1] from the expectation of the curve's displacement between [k -1, k]:

w * (x, t) = dt ∆t (x -E(ψ(x, k)|C k )), ∀t ∈ [k, k + 1], , (19)
where ∆t = t k+1 -t k and dt represents the time step associated to the dynamics discrete scheme. In this equation, E(•|C k ) denotes the expectation with respect to the path of C k up to time k (formally the natural filtration associated to the process C t ). This transportation velocity field is thus defined as a deterministic function computed from the realization of the curve (through its level set function ϕ t ) up to the previous frame instant k. As we will see it in the next section, this expectation will be empirically computed through a Monte Carlo approximation provided by our data assimilation filter. This approximation will be defined through a linear combination of the form:

w * (x, t) = dt ∆t x - 1 N N i=1 ω (i) ψ (i) (x, k) , ( 20
)
for t ∈ [k, k + 1]
, where ψ (i) corresponds to the auxiliary level set function associated to a given realization ϕ (i) of the implicit surface and ω (i) is a weight associated to the "quality" of this sample. Let us note that as w * (x, t) can only be used from instant t = 1, we have to define a transport component for t between 0 and 1. This initial transportation component is set to the velocity field estimated from the two first images through an optical-flow estimator.

To improve the relevance of this drift component, the velocity field is complemented with a local force derived from a segmentation process defined on color histograms [START_REF] Papadakis | A variational technique for time consistent tracking of curves and motion[END_REF]. This process 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 distributions. These two distributions are assumed to be estimated from the initial location of the target region. The segmentation is defined as the minimizer of the following energy function:

Ω F (ϕ, I)(x, t)dx, (21)
with

F (ϕ, I)(x, t) = d 2 B (ρ V x , ρ o )1 ϕ(x)<0 + d 2 B (ρ V x , ρ b )1 ϕ(x)≥0 , (22)
and where d 2 B is the squared Bhattacharya probability density distance measure and ρ V x is the probability density in the neighborhood V x of a pixel x. In our applications, the local neighborhood is defined as points such that y ∈ V x if |y -x| ≤ . By replacing the densities with intensity average, we retrieve the Chan and Vese functional proposed for image segmentation [START_REF] Chan | Active contours whithout edges[END_REF]. As we are using grayscale images with a pixel color range from 0 to 255, the Bhattacharya distance (which is a particular instance of the Hellinger distance between two probability distributions) is in our case defined as:

d 2 B (ρ 1 , ρ 2 ) = 1 2 255 0 ( ρ 1 (z) -ρ 2 (z)) 2 dz = 1 - 255 0 ρ 1 (z)ρ 2 (z)dz. ( 23
)
It is zero if the two densities are identical and equals to one for distributions with non-overlapping supports. The corresponding associated gradient operator with respect to the implicit surface in the sense of distributions is [START_REF] Chan | Active contours whithout edges[END_REF]:

∂ ϕ F = d 2 B (ρ V x , ρ o ) -d 2 B (ρ V x , ρ b ) δ(φ), (24)
where δ(•) is the Dirac function. An extension to all the level sets1 of this gradient is added to the drift previously defined from the vectorial function ψ [START_REF] Gordon | Novel approach to non-linear/non-gaussian bayesian state estimation[END_REF]. A timedependent linear combination of both components with proportions β(t) ∈ [0, 1] and 1 -β(t) reads then

w n = β(t)n T w * + (1 -β(t))∂ ϕ F (ϕ). ( 25
)
Let us remark, that we could perfectly add in the segmentation functional a term penalizing the curve length as it is usually the case in image processing:

σ τ Ω δ(ϕ)|∇ϕ(s)|ds. ( 26
)
Such a term would add a mean curvature component σ τ (ϕ)∇• ∇ϕ |∇ϕ| to the gradient (24) and we would then face a stochastic level set evolution equation of form ( 13), [START_REF] Dambreville | Tracking deformable objects with unscented kalman filtering and geometric active contours[END_REF] or [START_REF] Moral | Feynman-Kac Formulae Genealogical and Interacting Particle Systems with Applications[END_REF] depending on the balance between the noises diffusion as explained in the previous section.

The linear combination [START_REF] Merriman | Motion of multiple junctions: A level set approach[END_REF] enables to take into account the past trajectory of the curve and in the same time to apply corrections based on local color information. Within the interframe time range, the photometric component is likely to be especially helpful in the temporal vicinity of the second image, whereas the velocity component is more likely to be meaningful in the temporal vicinity of the first image. This remark allows us to define as a rule of thumb the proportion β of each component. In this study we choose to gradually change the proportion of each components according to:

β(t) = k + 1 -t, t ∈ [k, k + 1], (27)
so that in the one hand the photometric component is null at the beginning of the interframe time range and, in the other hand the past curve trajectories are not taken into account when the second image becomes available. This whole drift component together with the curve evolution law will allow us to draw sample of curves and to implement a stochastic image assimilation. To that end we will rely on Monte Carlo implementation of the filtering equations. We recall hereafter the principles governing such techniques.

Stochastic filtering, filtering distribution and Bayesian recursive filter

Stochastic filters aim at estimating the posterior probability distribution (called the filtering distribution) of a state variable trajectory (here the implicit surface ϕ 0:k ) starting from an initial state ϕ 0 ∈ R n up to the current time k given a complete sequence of measurements y 1:k = {y i , i = 1, . . . , k, y i ∈ R p } (that will be precised latter on for sake of clarity):

p(ϕ 0:k |y 1:k ). ( 28
)
At each discrete time instant k a measurement y k ∈ R p of the state is assumed to be available. These observations and the state variables are assumed to be link through a measurement equation:

y k = H(ϕ k ) + γ k , (29)
where γ k , the observation noise, is a white Gaussian noise with covariance matrix R, and operator H stands for the linear/nonlinear mapping from the state variable space to the observation space. The state variable at instant k is obtained by an integration of the level set evolution law

ϕ k = ϕ k-1 + k k-1 dϕ s ds, (30)
where the differential at a fixed grid point dϕ t (11) has a drift component defined through [START_REF] Merriman | Motion of multiple junctions: A level set approach[END_REF]. It is important to outline that as the level set evolution law is defined from an Ito integral, ϕ t is a Markov process. However, as the dynamics is nonlinear and involves in addition a multiplicative noise, the transition probability distribution p(ϕ k |ϕ k-1 ) is non Gaussian. Let us also point out that the numerical integration of ( 30) is performed with a time step δt that is usually about 10 to 100 times smaller than the latency δk between two subsequent measurements.

A recursive expression of the filtering distribution p(ϕ 0:k |y 1:k ) can be obtained from Bayes' theorem and the assumption that the measurements depend only on the current state (e.g. From Markoviannity of the dynamic system, we finally get the sought recursive expression:

p(ϕ 0:k |y 1:k ) = p(ϕ 0:k-1 |y 1:k-1 ) p(y k |ϕ k )p(ϕ k |ϕ k-1 ) p(y k |y 1:k-1 ) . (31)
For linear models, the filtering distribution is Gaussian and this recursive equation can be solved exactly through the well known Kalman filter, which provides a complete description of the pdf through explicit expressions of the two first moments.

In the case of nonlinear models, it is impossible to provide such an expression for the whole pdf. It is nevertheless possible to approximate the posterior distribution with a set of samples (called particles). This is presented in the next section.

Particle implementation of the nonlinear filtering

For nonlinear dynamics or nonlinear measurement models, a direct sampling from the filtering distribution is impossible since it would require the complete knowledge of the filtering distribution -which is in the general case non Gaussian -at a previous time.

Particle filtering techniques introduce a discrete approximation of the sought density as a sum of N weighted Diracs:

p(ϕ 0:k |y 1:k ) ≈ N i=1 ω (i) k δ ϕ 0:k (ϕ 0:k ) , (32)
centered on hypothesized locations of the state space sampled from a proposal distribution, π, (also called the importance distribution) approximating the true filtering distribution p(ϕ 0:k |y 1:k ). Each sample is then weighted by a coefficient ω

(i)
k , that accounts for the ratio between the two distributions for that sample. A great variety of importance functions can be chosen (with the only restriction that its support contains the filtering distribution one). Obviously the closer it is to the targeted filtering distribution the better the results. Under a weak hypothesis assuming a recursive factorization form of the importance function, the importance ratio can be recursively defined as:

ω (i) k ∝ ω (i) k-1 p(y k |ϕ (i) k )p(ϕ (i) k |ϕ (i) k-1 ) π(ϕ (i) k |ϕ (i) 0:k-1 , y 1:k ) . ( 33
)
By propagating the particles from time k -1 through the proposal density, and by weighting the sampled states with ω

(i)
k , a sampling of the filtering law is obtained. When the proposal distribution π is set to the dynamics π(ϕ

(i) k |ϕ (i) 0:k-1 , y 1:k ) = p(ϕ (i) k |ϕ (i) k-1 ), (34)
the weights updating rule [START_REF] Pons | Maintaining the point correspondence in the level set framework[END_REF] simplifies and implies only the data likelihood:

ω (i) k ∝ ω (i) k-1 p(y k |ϕ (i) k ). ( 35
)
This particular instance of the particle filter is called the Bootstrap filter or sequential importance resampling (SIR) filter [START_REF] Gordon | Novel approach to non-linear/non-gaussian bayesian state estimation[END_REF]. This is the kind of filter on which we will rely in this study. It is important to notice that for this filter, the direct sampling from the dynamics does not allow taking into account the current observation. The samples depend also only weakly on the past data through the filtering distribution estimated at the previous instant. This strategy in a general setting is thus completely "blind" to the measurements and requires a large number of particles to cover meaningful areas of the state space. This deficiency makes such a filter inoperant for high dimensional problem [START_REF] Snyder | Obstacles to high-dimensional particle filtering[END_REF]. To alleviate this problem it is very important to devise either an importance function that enables focusing on the most meaningful areas of the state space or to define a dynamics enabling to take into account the current measurements. In this study we choose the latter option, in defining the drift component from the curves trajectories on the previous interframe interval, but also from the current measurements. Let us note this model is not anymore a classical hidden Markov model as one of its component is driven by the data. It has been shown, however, that standard derivations can still be conducted with such models, leading to so-called conditional filters [START_REF] Arnaud | Partial linear Gaussian model for tracking in image sequences using sequential Monte Carlo methods[END_REF]. As a last ingredient, a resampling step of the particles is necessary to avoid the increase over time of the weight variance. This procedure discards particles with weak weights, and duplicates particles with high weights. For further elements concerning various aspects of particle filters (convergence properties, importance distribution choices, resampling strategies, etc.) interested readers may consult the textbooks [START_REF]The Oxford Handbook of Nonlinear Filtering[END_REF][START_REF] Moral | Feynman-Kac Formulae Genealogical and Interacting Particle Systems with Applications[END_REF][START_REF] Gordon | Sequential Monte Carlo methods in practice[END_REF].

Without entering into details, few words can be said on convergence. It has been shown by several authors that the root mean square error of the filtering distribution converges up to a constant as the inverse of the square root of the number of particles (C/ √ N ). The rate of convergence is hence quite slow and independent of the state space dimension. However, to reach a given accuracy on the mean square error the constant may depend on the state space dimension (see [START_REF] Crisan | Survey of convergence results on particle filering methods for practitioners[END_REF]). As a consequence the number of particles to consider may depend on the state space dimension. There is unfortunately to date no theoretical bound allowing to fix in practice the amount of particles. The number of particles considered is thus usually fixed on an experimental basis where a balance between the computational load of the method and a non diverging behavior of the filter is sought. For systems of large dimension implemented with only very few particles, due to evident computational limitations a residual error on the results will always remain. This error has obviously strong repercussions on the computation of all the moments of the filtering distribution.

Likelihood definition

In particle filtering procedures, the sought filtering distribution is expressed as a linear combination of samples [START_REF] Paragios | Geodesic active regions: a new framework to deal with frame partition problems in computer vision[END_REF]. The associated normalized weights, expressing the proportion of each sample in this discrete sum, depend on the data likelihood [START_REF] Pons | Maintaining the point correspondence in the level set framework[END_REF]. This is particularly true for the bootstrap filter, where the weights depends solely on this likelihood [START_REF] Sethian | Level set methods[END_REF]. Its definition is thus of crucial importance. A relevant choice of the likelihood enables penalizing efficiently curves that are too far from the expected results and to focus on meaningful areas of the state space.

In this work, we chose to define the likelihood from a distance between a color distribution of reference, ρ ref , and the color distribution, ρ (i) k , on the region delineated by the current curve sample. More precisely, for the Bhattacharyya distance, d B , [START_REF] Kimmel | Tracking level sets by level sets: a method for solving the shape from shading problem[END_REF], we define the data likelihood as:

p(y t |ϕ (i) t ) ∝ exp -λd 2 B (ρ ref ,ρ (i) k ) , ( 36
)
where λ is a parameter, that has been empirically set in practice to λ = 2. The reference color distribution can be defined in several ways. The simplest choice consists in specifying it from the empirical color histogram associated to an initial contour draw by the user. This is the procedure we employed. The likelihood and the weight computation can be finally summarized as follow:

-At time t = 0, a three-dimensional histogram is created, listing the colors contained within the curve.

-

At time t = k, an histogram h (i) k is computed for each particle ϕ (i) k (with i = 1, ..., N ). -The distances d (i)
B k between the initial histogram h 0 and the histograms h

(i) k are calculated -The likelihood v (i) k = e -λ(1-d (i) B k
) /Z and the weights ω

(i) k = ω (i) k-1 v (i) k /S, are finally computed with normalization factors Z = j v (j) k
and S = j ω (j) k respectively.

Parameters estimation

The stochastic dynamics proposed for the closed curve, depends on two diffusion parameters σ n and σ τ that weight the noises, respectively along the curve normal and tangent directions. In order to get the most efficient tracking as possible with only a very limited number of samples compared to the state space dimension, it is essential to consider diffusion coefficients that depend on time. As a matter of fact, the stronger the diffusion coefficients the wider the spread of the curve samples ensemble. This situation has to be engaged when the estimate or the data are of bad quality -e.g. low likelihood -, requiring thus to sustain an exploration of sufficiently large regions of the state space. At the opposite, when the current estimate is of good quality, one must stay in a close surrounding state space region and therefore set small diffusion values for the noises. Setting those coefficients from the data only may reveal quite involved, as it may require the introduction of some detection mechanisms. A simplest idea consists to fix them from the displacement field between two consecutive frames. This scheme is detailed is detailed hereafter. Let us note, u(x, k), the displacement undergoes by a curve's point from location x ∈ Ω at time k to its corresponding position at time k + 1. This displacement corresponds to a random path of the implicit function, ϕ, conditioned on the past observations. Assuming, this displacements field for the whole spatial domain constitutes a noisy version of the true displacement, w * , such that,

u(x, k) = w * (x, k) + σ n,k k+1 k ndB n,t + σ τ,k k+1 k n ⊥ dB τ,t , (37)
where the noises dB n,t and dB τ,t are directed along the normal and the tangent of the ϕ level-lines respectively. Besides, the diffusion coefficients are assumed constant on the time interval [t k , t k+1 ] and identical to those of the Ito diffusion [START_REF] Allaire | Structural optimization using sensitivity analysis and a level-set method[END_REF]. We are therefore making here the hypothesis that the noises associated with the level set displacement and the curve noises are collinear and have the same variances. Furthermore, we assume that the transport velocity field is such that

w * (x, t) = E(u(x, k)) dt ∆t , ∀t ∈ [k, k + 1].
The empirical covariances with respect to the filtering law of this observed displacement along the curve normal and tangent provide an estimation of the noise variances σ 2 n and σ 2 τ : and floods, which might afflict this area. Their analysis and forecasting are thus of the utmost interest for meteorologists and climatologists. As a first result, figures 1 and 3 plot four images of the evolution of the curve initialized by an ellipse on the first image of the infrared sequence (image #0 of figure 1). The second raw of both figures shows the trajectory along time of some of the curve's points. As can be observed from these two figures, the points trajectory recovered through the integration of the vectorial level function ψ attached to the curve's point location is nonlinear and strongly irregular at some places. The initial and the current curve location are plotted in dotted lines (respectively in red and green color) on the bottom row of this figure. On the top row, the red curve represents a curve sample curve. One example of results obtained for the whole sequence are plotted in figure 3. In this sequence the time latency between two satellite images is 15min.

σ 2 n,k = 1 N -1 N i=1   ω (i) C (i) u (i) (s, k) -w * (s, k) • n(s, t) 2 ds   , (38)
σ 2 τ,k = 1 N -1 N i=1   ω (i) C (i) u (i) (s, k) -w * (s, k) • n ⊥ (s, t) 2 ds   . (39)
To highlight the stability of the results with respect to the initial condition, we show on this example the results recovered for different runs. More precisely in this experiment, we run four times the algorithm for the very same initial condition. The results are plotted with different colors in figure 4 and5. Together with the four curves representing the zero level set of the empirical mean (with respect to the filtering distribution) we superimposed for the four curves an associated uncertainty band computed from the unbiased empirical variance of the implicit surface on the mean curve points (with respect to the filtering distribution): This variance defines the width of the uncertainty band. This band provides us a practical inspection tool enabling to visualize immediately the regions where the curve's samples exhibit some variability. This variability depends on the number of particles used to compute the empirical mean but also on the data. The spread is high in regions where the likelihood is poorly defined like in regions of low photometric gradient or far away from the target photometric distribution. A parallel situation occurs for motion estimation problems in low photometric gradient regions subject to high "aperture problems" where the uncertainty on the motion measurements is high (as on those regions the estimates highly depend on the prior distribution). Obviously, as deterministic procedures related to the Maximum a posteriori estimator (MAP) are used in practice, this uncertainty is not plainly revealed. The MAP may appear to be a very poor estimator for non linear possibly multimodal distribution. What should be in that case the result's representation? We believe a representation based on the two first moments is quite appropriate. However, for spaces of large dimension an empirical evaluation of those moments requires a huge amount of ensemble members (or realizations) that cannot be afforded computationally. A tradeoff must be found between a realistic number of particles and the results accuracy, with a strict limit provided by the filter divergence.

σ 2 (x) = 1 N -1 N i=1 ω i ϕ 2 i (x)δ(Eϕ i (x)).
To check if an amount of 100 particles was sufficiently high, we compared the results obtained for the four different runs with 100 and 200 particles respectively. These results are plotted figures 4 and 5. As can be observed, in agreement with the convergence results, the variations between the estimated curves is lower when the number of particles augments. The uncertainty bands associated to the curves are also superimposed. The different estimated curves remain globally close to each other although departures can be observed locally in regions associated with large uncertainty bands. These local differences depend on the data but also on the number of particles. For a larger number of particles the uncertainty band on those areas is larger. This seems to indicate that these uncertainties are mostly due to a locally poor likelihood information. As a matter of fact in such a situation all the curves proposed have locally very close probability and the spread increases with the number of particles in those regions. The distance between the curves however decreases clearly when the number of particles augments. We note also that all the curves keep the same topology in both cases.

To provide a quantitative criterion of the differences between the four different runs, we plotted on figure 5 several curves showing the evolution along time of the two first moments of the distance between the estimated curve together with the uncertainty band width square root (the variance is of the order of a squared distance for an implicit function taken in first approximation as a signed distance). The distance between the estimates is computed with a signed distance function applied on a given curve taken as a reference:

d(ϕ ref , ϕ i )(x) = |ϕ ref (x)| × δ(ϕ i (x)).
The empirical means and variances are then computed along the different estimated curves. As can be observed the mean and the variance of the distance between curves are lower for 200 particles. The accuracy is increased. However for 100 particles the mean distance is around or above 1 pixel for only two images (2.78 and 0.99 pixels respectively). The average distance for 100 particles is 0.51 pixels; for 200 particles it decreases to 0.35 pixels. Let us note the high mean and variance associated to the frame #42 is due to the disappearance of one of the curves one frame earlier on a region showing the evanescence of a convective cell. This causes a much higher distance between the curves in this region. For 200 particles all the curves are well in phase for this event and the distance is lower (see figure 6 for a quick comparison). For the two amounts of particles the filter does not show any error accumulation along time and demonstrates a stable behavior in term of the two first moment. The shape of the mean and variance of the uncertainty band width square root along the curves (i.e. the standard deviation of the level set function on the mean curves points) have close profiles. The variance is at some places slightly higher for 200 particles as explained previously. An amount of 100 particles constitutes a good tradeoff between accuracy and computational load. An amount of 200 particles increases only slightly the filter accuracy (the two first moment stability with respect to the initial condition).

As a second example, we show results obtained on a 71 frames Weather Watch Radar sequence in Fig. 8 with an interframe latency of 5min. With such sequences one aims at tracking areas associated with a high water density in order to detect the presence of rain and hail. The initial region is here set through a simple thresholding procedure of the areas with high water density. The successive locations of the curve tracked are plotted in figure 8. An uncertainty band around the estimate figure is also figured in gray color. This uncertainty band is given by the weighted local variance of the level set functions around the mean level set : the lower the variance, the narrower the uncertainty band. The tracking procedure proposed here allows us to track in a robust way the initial area of interest on a period of 350 minutes. Besides an evident interest for monitoring purpose, the system may be also of interest for short time forecasting applications. In order to

Number of particles: 200

Variance of distance between runs

Image index Number of pixels

Variance of standard deviation on the curve points (d) highlight this point we show on figure 9 and 10 for two different instants short time prevision at t + 25min. For this forecasting application the drift is fixed at the initial forecast instant and not changed during the whole prevision. There is thus no interaction with the data during this forecast. Examples of two previsions of 25 minutes are shown in figures 9 and 10. In those figures the yellow curve corresponds to the curve at the initial time of the forecast whereas the red curve is the forecasted curve. Although it is not used, we plotted also the radar image at the corresponding time to assess the quality of the prediction. In both cases, we can see that our drift component allows us to recover with a quite good agreement local changes of the region of interest. As a last example, we show the results that have been obtained on satellite images depicting the ice density of sea ice at the north pole. Those images are daily available. During summer it must be pointed out that due to the presence of an important sheet of water above the ice the measurements are completely unreliable. Besides temporal failures of the satellite communication system may generate large regions with no data. These failures may occur for several days. Despite such difficulties we run the method on several years of data in a complete autonomous way -without any hand tuning of parameters apart from the initial condition set the 09/09/1999 as a region with high ice density. A storyboard of the evolving curve during three years is shown in figures 11, 12 and 132 . As can be observed, the method succeeded to track reasonably well the target region along time. We see that during summer the curve tends to shrink near the north pole and then recover very quickly a meaningful area when reliable data becomes again available at the beginning of autumn. The uncertainty band corresponds also quite well to area of the curves for which the data does not correspond to the target ice distribution (through the photometric histogram). In figure 12 we focus on typical example of measurement failures. The behaviour of the method before, during and just after a satellite failure is plotted as in figures 12(g), 12(h), 12(i) and in figures 12(n), 12(o), 12(p), 12(q). This robust behaviour and the capacity to estimate the uncertainty associated to a given estimate is a salient feature of the method. Such abilities should be of interest in various applications in which one aims at tracking along time a region of interest from a noisy image sequence.

Conclusion

In this paper we have proposed a stochastic filtering framework that enables to track from an image sequence a closed curve represented as an implicit level set surface. This framework is built upon a stochastic representation the level set dynamics that relies on two different uncertainties on the curve motion, directed respectively along the curve normal and along the curve tangent. It also includes a transport vector field that combines an image-based force (related to local photometry) and a velocity induced by the past curve displacements and deformations. This deformation field enables the method to be robust to missing data along a long period of time. Besides, it has shown to provide a short time forecast model of good quality on Weather Watch radar data. Several examples of very different natures have demonstrated the generericity of the approach. In this work we focused on the tracking a single level set function, nevertheless tracking of segmentation maps through multiphase level sets [START_REF] Merriman | Motion of multiple junctions: A level set approach[END_REF][START_REF] Vese | A multiphase level set framework for image segmentation using the Mumford and Shah model[END_REF][START_REF] Zhao | A variational level set approach to multiphase motion[END_REF] could be probably envisaged within a similar framework. Extension to noises living in higher dimensional spaces could also be explored. As a by-product of this method, the transportation field could be informatively analysed or used to construct reduced order dynamics for interfaces separating complex phenomenon [START_REF] Artana | Strong and weak constraint variational assimilations for reduced order fluid flow modeling[END_REF].

Apendix

This section details the different terms involved in the expression of the evolution equation associated to the vectorial level set ψ.

Drift computation details

ψ is driven by the same velocity fields (5) as ϕ. For a fixed point y we have: In the same way as for ϕ, the differential of ψ(X t , t) at point X t = x reads:

dψ i t (y) = b i (
dψ i (X t ,
1 2 A i + F i + G i , (44)
with

A i = σ 2 n ∇ϕ T ∇

  p(y k |ϕ 0:k , y 1:k-1 ) = p(y k |ϕ k )): p(ϕ 0:k |y 1:k )p(y k |y 1:k-1 ) = p(ϕ 0:k , y k |y 1:k-1 ) = p(y k |ϕ 0:k , y 1:k-1 )p(ϕ 0:k |y 1:k-1 ) = p(y k |ϕ k )p(ϕ k |ϕ 0:k-1 , y 1:k-1 )p(ϕ 0:k-1 |y 1:k-1 ).

Fig. 1

 1 Fig. 1 Example of extracted convective cells from an infra-red satellite image sequenceimages #0 to #25. The second raw shows point trajectories associated to one sample curve.

Fig. 2

 2 Fig. 2 Example of extracted convective cells from an infra-red satellite image sequenceimages #33 to #45. The second raw shows point trajectories associated to one sample curve.

Fig. 3

 3 Fig. 3 Example of the tracking of moving convective cell, on a sequence of 81 satellite images of the infra-red channel, with 15 minutes between each image. The initial curve is initialized with ellipse on the first image

Fig. 4 Fig. 5

 45 Fig. 4 Example of four runs with 100 particles of the curve tracking filter on the infra-red atmospheric satellite images. All the curves have been identically initialized. The uncertainty band associated with the different curves are superimposed as transparent gray bands; the uncertainties are computed from the implicit surface variance on the zero level set of the filtering mean curve.

Fig. 6 Fig. 7

 67 Fig. 6 Evolution along time of the average and variance of (i) the distance between estimated curves obtained through several runs (curves plotted in blue color) and (ii) the uncertainty band width (in red color) for 100 and 200 particles respectively -the out of range value of the blue colored curve in graph (b) is 38.75.

Fig. 8

 8 Fig. 8 Example of tracking results on a Weather Watch Radar sequence. The initial curve is initialized around a region of high water density. The mean curve is plotted in red and the uncertainty band is displayed as a gray transparent band.

Fig. 9

 9 Fig. 9 Short time prevision of Weather Watch Radar region are shown each five minutes from t + 5 mn to t+25mn. The forecast curve at time t + k is plotted as a red curve whereas the curve at time t = 17 is displayed in yellow color. All the forecast curves depends only on the image data up tu time t.

Fig. 10

 10 Fig. 10 Short time prevision of Weather Watch Radar region are shown each five minutes from t + 5 mn to t+25mn. The forecast curve at time t + k is plotted as a red curve whereas the curve at time t = 28 is displayed in yellow color. All the forecast curves depends only on the image data up tu time t.

Fig. 11

 11 Fig. 11 Example of a region tracking on an ice density satellite sequence -year 1999 -2000.

Fig. 12

 12 Fig. 12 Example of a region tracking on an ice density satellite sequence -year 2000 -2001.

Fig. 13

 13 Fig. 13 Example of a region tracking on an ice density satellite sequence -year 2001 -2002.

  y, t)dt + f i (y, t)dB n,t + g i (y, t)dB τ,t .

	(40)

  2 ψ i t ∇ϕ + σ 2 τ (∇ϕ ⊥ ) T ∇ 2 ψ i t ∇ϕ ⊥ , F i = σ n |∇ϕ| ∇ϕ T ∇ 2 ψ i t ∇ϕ + ∇ 2 ϕ∇ψ i t -1 |∇ϕ| 2 (∇ψ i t ) T ∇ϕ∇ 2 ϕ∇ϕ , G i = σ τ |∇ϕ| ∇ϕ ⊥ T ∇ 2 ψ i t ∇ϕ ⊥ + ∇ 2 ϕ∇ψ i t -1 |∇ϕ| 2 (∇ψ i t ) T ∇ϕ ⊥ ∇ 2 ϕ∇ϕ .

								(45)
	The differential of ψ for a fixed point reads finally:
	dψ i t (x) = -(∇ψ i t ) T ∇ϕ |∇ϕ|	w * n dt -σ n (∇ψ i t ) T ∇ϕ |∇ϕ|	dB n,t -σ τ (∇ψ i t ) T ∇ϕ ⊥ |∇ϕ|	dB τ,t
	-	A i dt 2|∇ϕ| 2 +	σ n F i dt |∇ϕ|	+	σ τ G i dt |∇ϕ|	.	(46)

A molified version of the Dirac function[START_REF] Chan | Active contours whithout edges[END_REF][START_REF] Zhao | A variational level set approach to multiphase motion[END_REF] or simply δ = 1 can be used for this purpose.

A movie showing the method's results from 01/01/2005 to 18/11/2009 is also available as a supplementary material accompanying this paper.

Acknowledgments

We thank the CERSAT/IFREMER laboratory and Meteo-France for having provided us the ice density satellite image and the meteorological infra-red sequence. The authors acknowledge the ANR projects PREVASSEMBLE (ANR-08-COSI-012) and Geo-Fluids (ANR-09-SYSC-005) for their financial support.

For the particular case of the time interval between the two first images -where we do not have yet a first curve displacement estimation -, the values of these parameters are computed from the empirical variances of the initial motion field along the normal and tangent of the curve used to initialize our filter.

With those different ingredients, the whole method can be summarized by the algorithm (1). In the following section, we present the performances of our Algorithm 1 The Level-Set SIR filter, one iteration.

Require:

Set of level set samples at instant k -1: {ϕ

k-1 , i = 1, . . . , N } Diffusion coefficients values at instant k -1: σn(k -1) and στ (k -1) Drift at instant k -1: w * (x, k -1) for all x ∈ Ω Images I k and photometric histograms inside and outside the initial curve (ρo, ρ b ) Ensure:

Set of samples at time k: {ϕ

Experiments and results

This section reports several tracking results obtained with our approach. We aim here at highlighting the main abilities of the method. Note that in all these experiments, the number of particles is fixed as N = 100. As for the discrete schemes implemented, we choose a simple Euler Maruyama scheme for the temporal integration of the evolution laws associated to ϕ and ψ. For both evolution laws we rely on centered finite differences for the second order operators and on weighted ENO schemes [START_REF] Jiang | Efficient implementation of weighted eno schemes[END_REF][START_REF] Shu | Advanced Numerical Approximation of Nonlinear Hyperbolic Equations[END_REF] for the advection terms.

The first sequence of results concerns a Meteosat-MSG2 satellite image sequence of the infrared channel depicting the evolution of a convective cloud system. The analysis and tracking of such cloud systems is of major interest for aircraft and short time forecasters as they are known to be associated with hazardous meteorological events, such as strong wind drafts, lightnings, heavy rainfalls, hails or even tornadoes. Over tropical areas such as central Africa, convective cells produce most of the rain during the monsoon period, and are indirectly linked to droughts