É Debreuve

M Gastaud
email: gastaud@i3s.unice.fr

M Barlaud
email: barlaud@i3s.unice.fr

G Aubert
email: gaubert@math.unice.fr

Using the shape gradient for active contour segmentation: from the continuous to the discrete formulation

Keywords: Segmentation, shape optimization, shape derivative, active contour, discrete formulation, variational approach

A variational approach to image or video segmentation consists in defining an energy depending on local or global image characteristics, the minimum of which being reached for objects of interest. This study focuses on energies written as an integral on a domain of a function which can depend on this domain. The derivative of the energy with respect to the domain, the so-called shape derivative, is a function of a velocity field applied to the domain boundary. For a given, non-optimal domain, the velocity should be chosen such that the shape derivative is negative, thus indicating a way to deform the domain in order to decrease its energy. Minimizing the energy through an iterative deformation process is known as the active contour method. In the continuous framework, setting the velocity to the opposite of the gradient associated with the L 2 inner product is a common practice. In this paper, it is noted that the negativity of the shape derivative is not preserved, in general, by the discretization of this velocity required by implementation. In order to guarantee that the negativity condition holds in the discrete framework, it is proposed to choose the velocity as a linear combination of pre-defined velocities. This approach also gives more flexibility to the active contour process by allowing to introduce some a priori knowledge about the optimal domain. Some experimental results illustrate the differences between the classical and the proposed approach.

Introduction

The shape gradient approach is presented in Sections 1 and 2 and Appendices A and B. Apart from presenting the context, these sections emphasize the fact that 1.1 A variational approach to segmentation

Image or video segmentation can be performed with the following variational approach: the problem is formulated as the minimization of an energy depending on assumed characteristics of objects of interest [START_REF] Caselles | Geodesic active contours[END_REF][START_REF] Chesnaud | Statistical region snake-based segmentation adapted to different physical noise models[END_REF][START_REF] Cremers | Variational space-time motion segmentation[END_REF][START_REF] Jehan-Besson | DREAM 2 S: Deformable regions driven by an eulerian accurate minimization method for image and video segmentation[END_REF][START_REF] Li | Dual front evolution model and its application in medical imaging[END_REF][START_REF] Li | Local or global minima: flexible dual front active contours[END_REF][START_REF] Mumford | Optimal approximations by piecewise smooth functions and associated variational problems[END_REF]. For simplicity, it is supposed that there is a unique object. Typically, the energy is a sum of a domain integral and a boundary integral [START_REF] Gastaud | Combining shape prior and statistical features for active contour segmentation[END_REF][START_REF] Mumford | Optimal approximations by piecewise smooth functions and associated variational problems[END_REF][START_REF] Paragios | Geodesic active regions for motion estimation and tracking[END_REF]]

E(Γ) = Ω φ f (Γ, x) dx + Γ ϕ f (s) ds ( 1
)
where Ω is an open set of R2 , Γ is the oriented boundary ∂Ω of Ω, s is the arc-length parameterization of Γ, and f is the image or video to be segmented1 . Energy (1) is designed to have a unique global minimum 2 at Ω ⋆ , the domain of the object of interest. For an ideal object, if φ is equal to zero, φ is sometimes referred to as the descriptor of the object. For instance, if φ is equal to

φ(Γ, x) = (f (x) -µ(Γ)) 2 (2)
where µ(Γ) is the average value of f in Ω, then φ is equal to zero on Ω if and only if f is constant on Ω. Therefore, φ is a descriptor of objects of constant intensity. Similarly, ϕ is the descriptor of the object boundary.

If the object background can also be characterized, then the energy can be symmetrized as follows

E(Γ) = Ω φ(Γ, x) dx + Γ ϕ(s) ds + Ω c φ c (Γ c , x) dx (3)
where Ω c is the complement D\ Ω of Ω in image or frame domain D. This combination of an integral on Ω and an integral on Ω c is called region competition [START_REF] Cremers | Motion competition: Variational integration of motion segmentation and shape regularization[END_REF][START_REF] Debreuve | Space time segmentation using level set active contours applied to myocardial gated SPECT[END_REF][START_REF] Yezzi | A statistical approach to snakes for bimodal and trimodal imagery[END_REF]. Nevertheless, in the following developments, the last integral of (3) will be discarded since its handling is similar to the one of the first integral.

Active contour

A possible approach to minimize energy (1) is to set an initial contour and to deform it iteratively in such a way that its energy decreases from one iteration to the next. The contour eventually converges toward a (possibly local) minimizer. This process is known as active contour [START_REF] Caselles | A geometric model for active contours[END_REF][START_REF] Caselles | Geodesic active contours[END_REF][START_REF] Chan | Active contours without edges[END_REF][START_REF] Kass | Snakes: Active contour models[END_REF]. When φ does not depend on Γ, the minimization procedure can be based on a local strategy involving energy evaluations only [START_REF] Germain | Optimal snakebased segmentation of a random luminance target on a spatially disjoint background[END_REF]. Otherwise, the influences of local deformations are linked together and the contour deformation should be performed at once. The appropriate deformation can be derived from the derivative of the energy with respect to the contour [START_REF] Caselles | A geometric model for active contours[END_REF][START_REF] Caselles | Geodesic active contours[END_REF][START_REF] Kass | Snakes: Active contour models[END_REF].

Shape derivative and evolution equation

The derivative of energy [START_REF] Arulampalam | A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking[END_REF] with respect to Γ can be obtained by a calculus of variations [START_REF] Aubert | Image segmentation using active contours: Calculus of variations or shape gradients?[END_REF][START_REF] Caselles | Geodesic active contours[END_REF]. However, if descriptors φ or ϕ depends on Γ, this can be complex. Some studies on shape gradients [START_REF] Delfour | Shapes and geometries: Analysis, differential calculus and optimization[END_REF][START_REF] Jehan-Besson | DREAM 2 S: Deformable regions driven by an eulerian accurate minimization method for image and video segmentation[END_REF][START_REF] Schnörr | Computation of discontinuous optical flow by domain decomposition and shape optimization[END_REF][START_REF] Sokolowski | Introduction to shape optimization: Shape sensitivity analysis[END_REF] offer a convenient and general basis for this differentiation.

The shape derivative of (1) is a function dE(Γ, V ) of Γ and a velocity V defined on Ω but applied to Γ only. For a given Γ, the velocity should be chosen such that the shape derivative is negative, thus indicating a way to deform Γ in order to decrease its energy. As noted by [START_REF] Charpiat | Designing spatially coherent minimizing flows for variational problems based on active contours[END_REF], the straightforward choice is to take the opposite of the gradient of (1) associated with the L 2 inner product on Γ. However, one might want to use other descent directions, for example, to improve the convergence rate [START_REF] Hintermüller | A second order shape optimization approach for image segmentation[END_REF], to increase spatial coherence and avoid as much as possible to converge to irrelevant local minima [START_REF] Charpiat | Designing spatially coherent minimizing flows for variational problems based on active contours[END_REF], to simplify implementation [START_REF] Overgaard | An analysis of variational alignment of curves in images[END_REF][START_REF] Solem | A geometric formulation of gradient descent for variational problems with moving surfaces[END_REF], or to respect the principles of an underlying physical model or for improved stability and convergence rate [START_REF] Doǧan | Finite element methods for shape optimization and applications[END_REF]. These alternatives may result from (i) designing other inner products or from (ii) directly designing descent directions V verifying [START_REF] Overgaard | An analysis of variational alignment of curves in images[END_REF][START_REF] Solem | A geometric formulation of gradient descent for variational problems with moving surfaces[END_REF] dE(Γ, V ) ≤ 0 .

In any case, the active contour evolution equation has the form ∂Γ ∂τ = V (τ ) .

When implementing this evolution equation, the negativity of the shape derivative is not preserved, in gen-eral, by the discretization of velocity V . This paper proposes to select other descent directions allowing to solve this problem in the context of (ii). Moreover, the proposed approach offers the possibility to introduce some a priori knowledge about the optimal contour. Section 2 reminds some results about the shape derivative. Section 3 describes the above mentioned discretization error corresponding to the so-called direct approach. It proposes a solution, called the constrained approach, based on choosing V as a linear combination of pre-defined velocities. Section 4 provides additional remarks about the constrained approach. Section 5 presents some experimental results to illustrate the differences between the direct approach and the constrained approach, and to illustrate how the constrained approach allows to introduce some a priori knowledge about the solution of (5).

Notations

In the following, D is a subset of R 2 and the image to segment, f , is a function from D to R m . Domain Ω is an open set of D. Γ is the oriented boundary ∂Ω of Ω, and s is the arc-length parameterization of Γ. For convenience, the notation a(s) refers to a(Γ(s)). Samples on Γ are denoted by Γ

i , i ∈ [1, n]. Arc-length s i is such that Γ(s i ) = Γ i .
Then, both notations a(s i ) and a(Γ i ) are equivalent. Note that s 1 is equal to 0 and s n is equal to L -(s n+1 -s n ) where L is the length of Γ and Γ(s n+1 ) = Γ(s 1 ). The contour segment between Γ i and Γ i+1 (a line segment, a spline segment. . . depending on the contour representation) is denoted by γ i .

The L 2 inner product on the space of velocities will be denoted

U, V = Γ U (s) • V (s) ds ( 6
)
where • is the dot product.

2 Expression of the shape derivative

General expression

Let us consider the typical energy

E(Γ) = Ω φ(Γ, x) dx + Γ ϕ(s) ds . (7)
Energy ( 7) is a function of a contour and the set of contours is not a vector space. As a consequence, the derivative of ( 7) with respect to Γ cannot be expressed in the usual way. Let Ω(τ ), τ ≥ 0, be a family of domains such that Ω(τ = 0) = Ω. Then, when τ increases, Γ(τ ) can be considered as a deforming interface in a medium characterized by φ and ϕ. Hence, a development inspired by works in continuum mechanics [START_REF] Haug | Methods of engineering mathematics[END_REF] can be made in order to determine the derivative of [START_REF] Charpiat | Designing spatially coherent minimizing flows for variational problems based on active contours[END_REF] with respect to τ at τ equal to zero [START_REF] Debreuve | Space time segmentation using level set active contours applied to myocardial gated SPECT[END_REF].

The study of such energies and their variations was further developed in the framework of shape optimization [START_REF] Delfour | Shapes and geometries: Analysis, differential calculus and optimization[END_REF][START_REF] Hintermüller | A second order shape optimization approach for image segmentation[END_REF][START_REF] Jehan-Besson | DREAM 2 S: Deformable regions driven by an eulerian accurate minimization method for image and video segmentation[END_REF][START_REF] Schnörr | Computation of discontinuous optical flow by domain decomposition and shape optimization[END_REF][START_REF] Sokolowski | Introduction to shape optimization: Shape sensitivity analysis[END_REF]. In this context, the following expression is known as the shape derivative of ( 7)

dE(Γ, V ) = Ω ∂φ(Γ(τ ), x) ∂τ τ =0 dx - Γ φ(Γ, s) - ∂ϕ(s) ∂N + ϕ(s) κ(s) N (s) • V (s) ds (8)
where vector field V is by definition a velocity defined on Ω, N is the inward unit normal of Γ, and κ is the curvature of Γ (some details are given in Appendices A and B).

Expression with only a boundary integral

Under weak assumptions, the shape derivative of the domain integral in [START_REF] Charpiat | Designing spatially coherent minimizing flows for variational problems based on active contours[END_REF] has an equivalent expression in the form of a boundary integral [START_REF] Sokolowski | Introduction to shape optimization: Shape sensitivity analysis[END_REF][START_REF] Solem | A geometric formulation of gradient descent for variational problems with moving surfaces[END_REF]. As is clear from Section 3.1.1, such an expression is convenient in the active contour framework. Indeed, the evolution equation of an active contour can be easily deduced from a boundary integral-only expression.

In the following, it is assumed that the shape derivative of [START_REF] Charpiat | Designing spatially coherent minimizing flows for variational problems based on active contours[END_REF] has been rewritten into a boundary integral

dE(Γ, V ) = - Γ Ψ(Γ, s) N (s) • V (s) ds, (9)
either because one of the two conditions (87) or (104) applies (see Appendix C), or as the result of another development.

3 From the continuous to the discrete framework

In this section, the classical, or direct, approach to using the shape derivative in the active contour framework is presented. It corresponds to choosing the gradient associated with the L 2 inner product as the descent direction (among the velocities that ensure the negativity of the shape derivative) and to discretizing it. It will be noted that this direct approach implies a discretization error possibly responsible for the loss of the negativity condition. A constrained approach will be proposed to solve this problem. The constraint is specified in terms of a linear combination of predefined velocities. It will be shown that the negativity condition still holds after discretization when choosing the pre-defined velocities in accordance with the active contour representation.

Direct approach

3.1.1 Negativity of the shape derivative Shape derivative ( 9) is a function of a velocity field V . Since energy [START_REF] Charpiat | Designing spatially coherent minimizing flows for variational problems based on active contours[END_REF] must be minimized, it is necessary to choose V such that ( 9) is negative. In other words, the following inequality must be solved for V

dE(Γ, V ) ≤ 0 . (10)
Interpreting ( 9) as the L 2 inner product on the space of velocities, the velocity

G = -Ψ N (11)
can be identified with the gradient associated with this inner product. It is called the shape gradient of (7) [START_REF] Charpiat | Designing spatially coherent minimizing flows for variational problems based on active contours[END_REF][START_REF] Delfour | Shapes and geometries: Analysis, differential calculus and optimization[END_REF][START_REF] Hintermüller | A second order shape optimization approach for image segmentation[END_REF][START_REF] Sokolowski | Introduction to shape optimization: Shape sensitivity analysis[END_REF]. Then, taking a steepest descent approach, one can think of choosing [START_REF] Debreuve | Space time segmentation using level set active contours applied to myocardial gated SPECT[END_REF][START_REF] Jehan-Besson | DREAM 2 S: Deformable regions driven by an eulerian accurate minimization method for image and video segmentation[END_REF] V

(s) = -G(Γ, s) (12)
in the following evolution equation

∂Γ ∂τ = V (τ ) . (13)

Minimization in the continuous framework

The implementation of evolution equation ( 13) can be based on a finite difference approximation of the derivative with respect to τ verifying the CFL-condition. Instead, a line search strategy can be used [START_REF] Hintermüller | A second order shape optimization approach for image segmentation[END_REF] Γ 0

Γ +1 = Γ + α V ( 14
)
where Γ 0 is an initial contour, superscript +1 represents the next element of a sequence3 , α is a positive constant, and V is given by [START_REF] Debreuve | Active contour segmentation with a parametric shape prior: link with the shape gradient[END_REF]. Note that V is defined on Ω (see Appendix B) but given on Γ only. However, as is clear from [START_REF] Delingette | Shape and topology constraints on parametric active contours[END_REF], it does not need to be known anywhere else. The optimal value for α can be computed as follows

α = arg min α≥0 E(Γ +1 (α)) . (15)

Discretization and induced velocity

In practice, active contour Γ is sampled. For instance, it can be represented by a polygon {Γ i , i ∈ [1, n]} without self-intersection [START_REF] Germain | Optimal snakebased segmentation of a random luminance target on a spatially disjoint background[END_REF]. The corresponding discrete version of evolution equation ( 14) is

Γ 0 i Γ +1 i = Γ i + α V (Γ i ) ( 16
)
where {Γ 0 i } is an initial polygon. Note that ( 16) does not make use of V along the edges γ i of the polygon. Instead, it implicitly defines a velocity Ṽ , called induced velocity, transforming edges γ i into edges γ +1

i . However, it is unlikely that such a transformation corresponds to the chosen velocity4 (12) (see Fig. 1a). As a consequence, the negativity of shape derivative (9) at Ṽ is not guaranteed. In other words, the discrete evolution equation ( 16) might not generate a minimizing sequence of contours. This problem is not specific to the polygonal representation. It also arises if the contour is represented by a smooth curve, say, a spline, since it is due to the sampling of V (see Fig. 1b). Generally speaking, only these contour samples where V is computed are translated correctly. The contour segments in-between the translated samples are defined a posteriori by the chosen interpolation model (a polygon, a uniform cubic B-spline. . . ).

However, as pointed out by a reviewer, if the edge length is small enough, dE(Γ, Ṽ ) is still negative and this direct approach is valid (see Appendix D). The proposed edge length upper bound is global and depends on maximal variations of Ψ over Γ. It could also be adapted locally to be low in portions of Γ where Ψ varies a lot and higher where Ψ varies slowly [START_REF] Taton | Deformable model with non-Euclidean metrics[END_REF].

Finally, note that (16) might converge too early since the condition {V (Γ i ) = 0, i ∈ [1, n]} is less restrictive than {V (s) = 0, s}. But again, the smaller the edge length, the less critical. Nevertheless, instead of a condition on the edge length, one can wonder if there is a way to choose V such that Ṽ is equal to V .

Constrained approach

Negativity of the shape gradient

In order to guarantee that, after discretization, the induced velocity Ṽ matches the original velocity V , we propose to restrict, beforehand in the continuous domain, the domain transformations5 to a linear combination of a set of pre-defined transformations

T (τ ) = i β i T i (τ ), β i ∈ R . ( 17
)
The differentiation of ( 17) with respect to τ leads to

V = i β i V i ( 18
)
where V i is a so-called pre-defined velocity. In this context, velocity V will be given by a choice of β i 's that satisfies the negativity of shape derivative [START_REF] Cremers | Motion competition: Variational integration of motion segmentation and shape regularization[END_REF]. The shape derivative can be rewritten as

dE(Γ, V ) = dE(Γ, i β i V i ) ( 19
) = i β i dE(Γ, V i ) (20) = β • d E(Γ, V i ) ( 21
)
where • is the dot product, β is the vector of components β i , and d E(Γ, V i ) is the vector of components dE(Γ, V i ). Taking a steepest descent approach, β should be such that dE(Γ, V ) is as negative as possible.

The Cauchy-Schwarz inequality implies that

|β • d E(Γ, V i )| ≤ |β| |d E(Γ, V i )| (22)
and, for a given norm of β, dE(Γ, V ) is as large as possible in absolute value when β and d E(Γ, V i ) are linearly dependent. Therefore, the β i 's are taken equal to

β i = -α dE(Γ, V i ) ( 23
)
where α is a positive constant. Section 3.2.3 presents a general purpose definition of pre-defined velocities and an application-driven variation of it.

Minimization in the continuous framework and discretization

In the continuous framework, a possible minimizing sequence is

Γ 0 Γ +1 = Γ -α i dE(Γ, V i ) V i ( 24
)
where the optimal value for α can be computed as follows α = arg min

α≥0 E(Γ +1 (α)) . ( 25
)
With the polygonal representation mentioned in Section 3.1.3, the discretization of ( 24) leads to the following evolution equation

Γ 0 i Γ +1 i = Γ i -α j dE(Γ, V j ) V j (Γ i ) . ( 26
)
However, for the pre-defined velocities proposed in Section 3.2.3, the sum over the pre-defined velocities reduces to a single term since at Γ i , only V i is not equal to zero

Γ 0 i Γ +1 i = Γ i -α dE(Γ, V i ) V i (Γ i ) . (27)

Examples of pre-defined velocities for the polygonal representation

According to the remarks of Section 3.1.3, it can be deduced that the pre-defined velocities must be consistent with the contour representation: the representation must be preserved when the contour is deformed by a linear combination of the V i 's. Moreover, ideally, the V i 's should be able to generate any velocity V by linear combination (at least, if there is no a priori knowledge about the optimal contour). There is no such basis of velocities. However, the V i 's must allow to represent a reasonable variety of velocities. If the contour is represented by a polygon with n vertices, the following definitions can be considered.

(i) At vertex Γ i , a pseudo-normal N (Γ i ) is de- fined [14, 27] and V i is the velocity collinear to N (Γ i ) at Γ i and transforming (Γ i-1 , Γ i , Γ i+1 ) into (Γ i-1 , (Γ i + N (Γ i )), Γ i+1 ) (see Fig. 2). In other words, V i is a vector field with support [Γ i-1 , Γ i+1 ], linear from zero (at Γ i-1 ) to N (Γ i ) (at Γ i )
, and linear again back to zero (at Γ i+1 ). This definition involves no a priori knowledge (see Section 5.1 for an example of application).

(ii) The previous definition can be modified to introduce some a priori knowledge. For instance, in tracking, the approximate motion of the tracked object might be known. In particular, a joint segmentation and motion computation method [START_REF] Cremers | Variational space-time motion segmentation[END_REF] certainly requires computation of the motion for a given, fixed segmentation (and vice-versa). Therefore, if m i is the estimated motion of Γ i , then V i can be defined similarly to (i) by replacing N (Γ i ) with m i or m i /|m i | (see Section 5.2 for an example of application).

Coherence between the continuous and discrete evolutions

In order to establish that ( 24) and ( 27) lead to identical evolutions, it must be shown that the deformations induced by ( 27) on contour segments γ i are equal to the deformations of the continuous evolution [START_REF] Kass | Snakes: Active contour models[END_REF]. This would also prove that the negativity of the shape derivative is preserved in the discrete framework. Indeed, if both evolutions are identical, the induced velocity Ṽ matches the constrained velocity V .

Polygonal representation Let x be a point on edge γ i (see Fig. 3). There exists t ∈ [0, 1] such that

x = γ i (t) = (1 -t) Γ i + t Γ i+1 . ( 28
)
Let V i be a pre-defined velocity as defined in Sec-

V (t) = β i V i (t) + β i+1 V i+1 (t) . ( 29
)
Therefore, in the continuous framework, point γ i (t) is translated to

γ i (t)+V (t) = (1-t) Γ i +β i V i (t)+t Γ i+1 +β i+1 V i+1 (t) . (30)
The pre-defined velocities are such that

V i (t) = (1 -t) V i (Γ i ) V i+1 (t) = t V i+1 (Γ i+1 ) . (31)
Combining ( 30) and ( 31) leads to

γ i (t) + V (t) = (1 -t) (Γ i + β i V i (Γ i )) +t (Γ i+1 + β i+1 V i+1 (Γ i+1 )) (32) = (1 -t) Γ +1 i + t Γ +1 i+1 , (33)
meaning that the translation of point γ i (t) in the continuous framework (left-hand side) is a point on line segment [Γ +1 i , Γ +1 i+1 ] obtained by joining the discrete translations of Γ i and Γ i+1 (right-hand side). Moreover, when γ i (t) describes edge γ i , then the translation of γ i (t) describes the whole line segment. As a conclusion, the continuous equation ( 24) and the discrete equation ( 27) lead to identical evolutions.

Uniform cubic B-spline representation

If the contour is represented by a uniform cubic B-spline with n segments [START_REF] Brigger | B-spline snakes: A flexible tool for parametric contour detection[END_REF][START_REF] Jacob | A unifying approach and interface for spline-based snakes[END_REF][START_REF] Precioso | Robust real-time segmentation of images and videos using a smooth-spline snake-based algorithm[END_REF][START_REF] Unser | B-spline signal processing: Part I -Theory[END_REF] (see Appendix E for the notations and properties used here), then one can think of basing the pre-defined velocities on blending function B (see Fig. 9). For example,

V i (t) = B i i (t) N (Γ i ), t ∈ [-n/2, n/2] (34)
where B i i represents blending function B shifted to be "centered" on Γ i and "seen" from spline segment γ i . More specifically,   

B i i (Γ i ) = B(0) B i i (Γ i-1 ) = B i i (Γ i+1 ) = B(1) B i i (Γ i-2 ) = B i i (Γ i+2 ) = B(2) (35)
and Γ j corresponds to t = j -i. Therefore, velocity V on spline segment γ i parameterized as described above is equal to

V (t) = j β j V j (t) (36) = i+2 j=i-1 β j V j (t) (37) = i+2 j=i-1 β j B i j (t) N (Γ j ) . (38)
Moreover, γ i is defined by

γ i (t) = i+2 j=i-1 q j B i j (t), t ∈ [0, 1] . (39)
Therefore, after evolution, point γ i (t) is translated to

γ +1 i (t) = γ i (t) + V (t) (40) = i+2 j=i-1 (q j + β j N (Γ j )) q +1 j B i j (t) . (41)
Equation ( 41) represents a spline uniquely defined by the new control points q +1 i . It must be checked whether the sampling points Γ +1 i this spline interpolates are also obtained by the discrete evolution equation. Since [START_REF] Yezzi | A statistical approach to snakes for bimodal and trimodal imagery[END_REF] represents a closed, uniform cubic B-spline, the following property holds

Γ +1 i = 1 6 (q +1 i-1 + 4q +1 i + q +1 i+1 ) (42) = 1 6 [(q i-1 + 4q i + q i+1 ) + β i-1 N (Γ i-1 ) +4β i N (Γ i ) + β i+1 N (Γ i+1 )] (43) = Γ i + β i-1 B i i-1 (0)N (Γ i-1 ) + β i B i i (0)N (Γ i ) +β i+1 B i i+1 (0)N (Γ i+1 ) (44) = Γ i + i+1 j=i-1 β j B i j (0) N (Γ j ) (45) Noting that B i i+2 (0) is equal to zero, (45) is equiva- lent to Γ +1 i = Γ i + i+2 j=i-1 β j B i j (0) N (Γ j ) (46) = Γ i + V (0) (47)
which is the discrete evolution equation [START_REF] Li | Local or global minima: flexible dual front active contours[END_REF]. Therefore, this spline representation associated with appropriate pre-defined velocities guarantees that the discrete evolution matches the continuous evolution [START_REF] Kass | Snakes: Active contour models[END_REF].

4 Some remarks about the predefined velocities

Normalization

The weight of pre-defined velocity V i is equal to -α dE(Γ, V i ). The weight of k V i is k times larger. Therefore, a pre-defined velocity can be made artificially preponderant in evolution process [START_REF] Li | Local or global minima: flexible dual front active contours[END_REF]. As a consequence, it seems appropriate to normalize the predefined velocities (unless otherwise imposed by some a priori). Actually, they are defined up to a multiplicative constant. Their norm does not need to be equal to one. It should only be the same for all V i 's. The pre-defined velocities proposed in Section 3.2.3 are such that

|V i | 2 = V i , V i = 2l 3 ( 48
)
where l is the edge length.

Interpretation as a projection

Individual projections and redundancy

Weights ( 23) can be interpreted as individual projections of the opposite of gradient G (see Eq. ( 11)) onto the V i 's

β i = -G, V i , (49)
up to a multiplicative constant. One can check that the pre-defined velocities proposed in Section 3.2.3 do not form an orthogonal set relative to the L 2 inner product

V i , V j =    2l 3 if i = j l 6 N (Γ i ) • N (Γ j ) if |i -j| = 1 0 otherwise (50)
where

N (Γ i ) • N (Γ j
) is equal to the cosine of the angle between N (Γ i ) and N (Γ j ). In other words, these predefined velocities are redundant.

Let us consider three successive vertices Γ i-1 , Γ i , and Γ i+1 and their respective pre-defined velocities. Suppose that G is equal to V i . If the pre-defined velocities did form an orthogonal set, the velocity V applied to Γ would be proportional to -V i , as expected. Instead, the pre-defined velocity weights are equal to (up to a multiplicative constant) -N (Γ i-1 ) • N (Γ i ), -4, and -N (Γ i ) • N (Γ i+1 ), respectively. It means that all three pre-defined velocities get weights different from zero (in general) with a higher weight for V i (in absolute value). If the dot products are positive (i.e., if the (pseudo) curvature is not too high in absolute value), this can be interpreted as a smoothing of the response obtained in the case of an orthogonal set of pre-defined velocities, which could, as suggested in [START_REF] Charpiat | Designing spatially coherent minimizing flows for variational problems based on active contours[END_REF], increase the spatial coherence of the active contour evolution. Unfortunately, if the dot products are negative, V i-1 and V i+1 get positive weights, leading to an undesirable evolution where Γ i is translated in a way that decreases the energy and the two other vertices are translated in a way that (although slightly) tends to increase the energy.

As suggested by a reviewer, the redundancy can be suppressed by generating mutually orthogonal predefined velocities V ⊥ i using the Schmidt orthonormalization process, or G could be projected onto the space of linear combinations of the pre-defined velocities instead of being projected individually onto each predefined velocity (see Section 4.2.2).

Alternative linear combination

Instead of computing weights β i as individual projections (49), they could be computed such that V is the L 2 projection of -G onto the space of linear combinations of the V i 's

β = arg min β -G - i β i V i 2 (51) = arg min β i β i V i 2 + 2 i β i G, V i (52)
= arg min

β i β 2 i |V i | 2 + 2 i<j β i β j V i , V j +2 i β i G, V i (53)
where β is the vector of components β i .

If the V i 's form an orthonormal set, then (53) is equivalent to

β = arg min β |β| 2 + 2 β • η ( 54
)
where η is the vector of component G, V i . For a given norm of β, the Cauchy-Schwarz inequality leads to the result based on individual projections β i = -η i , up to a multiplicative constant.

For the pre-defined velocities proposed in Section 3.2.3, if G is equal to V i , then the projection of -G is necessarily such that β j is equal to zero if j ∈ {i -1, i, i + 1}. Let us fix the norm of β to one. Using the fact that

β 2 i-1 + β 2 i + β 2
i+1 is equal to one, one can check that the projection is unique and given by

β i = -1 β j = 0, j = i , (55)
thus solving the problem of redundancy mentioned in Section 4.2.1. The constrained approach using the projectionbased weights (53) was implemented and compared with the constrained approach using the weights based on individual projections (49) (see Section 5.1).

Parametric approach

Presentation

Another way to avoid the problem mentioned in Section 3.1.3 is to parameterize the active contour and to rewrite energy [START_REF] Charpiat | Designing spatially coherent minimizing flows for variational problems based on active contours[END_REF] as a function of these parameters [START_REF] Jacob | Efficient energies and algorithms for parametric snakes[END_REF][START_REF] Marmorat | Private communication on shape gradient in the context of a parametric active contour[END_REF][START_REF] Unal | Informationtheoretic active polygons for unsupervised texture segmentation[END_REF]. The problem of minimizing [START_REF] Charpiat | Designing spatially coherent minimizing flows for variational problems based on active contours[END_REF] is then restricted to the set of domains whose boundary can be described by such parameters. It can be solved in R n where n is the number of parameters. If φ in [START_REF] Charpiat | Designing spatially coherent minimizing flows for variational problems based on active contours[END_REF] does not depend on Γ, then the partial derivatives of [START_REF] Charpiat | Designing spatially coherent minimizing flows for variational problems based on active contours[END_REF] with respect to each parameter can be obtained by a calculus of variations [START_REF] Jacob | Efficient energies and algorithms for parametric snakes[END_REF]. However, it can be more complex if φ does depend on Γ, unless noticing that the partial derivatives can be expressed as shape derivatives [START_REF] Debreuve | Active contour segmentation with a parametric shape prior: link with the shape gradient[END_REF][START_REF] Marmorat | Private communication on shape gradient in the context of a parametric active contour[END_REF]. Let Γ be a curve described by a set of parameters p = {p i , i ∈ [1, n]}, e.g., a spline. Energy [START_REF] Charpiat | Designing spatially coherent minimizing flows for variational problems based on active contours[END_REF] can be rewritten as a function of p

E(p) = Ω(p) φ(Γ(p), x) dx + Γ(p) ϕ(s) ds . (56)
It can be shown that the gradient of (56) is equal to

∇E(p) = i dE Γ, ∂Γ ∂p i e i ( 57
)
where ∂Γ/∂p i is a so-called admissible velocity and e i is the i th element of the canonical basis of R n . Energy (56) can be minimized using the following procedure

p 0 p +1 = p -α ∇E(p) (58)
where the optimal value for α can be computed as follows α = arg min α≥0

E(p +1 (α)) . (59)
Section 4.3.2 shows that the constrained approach ( 24)/( 26) is equivalent to the parametric approach (58) in certain circumstances.

Link with the constrained approach

Polygonal representation If the contour is represented by a polygon with n edges (see Section 3.2.3), then the parameters involved in the parametric approach are the coordinates (a i , b i ) of each vertex Γ i . Referring to the notations of Section 4.3.1, parameters p are

p = {p j , j ∈ [1, 2n]} = {a 1 , b 1 , a 2 , b 2 . . . , a n , b n } where (a i , b i ) = Γ i .
(60) Therefore, there are 2n admissible velocities such that

       ∂Γ ∂a i (Γ i ) = e 1 ∂Γ ∂b i (Γ i ) = e 2 , i ∈ [1, n] (61)
where (e 1 , e 2 ) is the canonical basis of R 2 . Writing procedure (58) in terms of vertices leads to

Γ +1 i = Γ i -α dE Γ, ∂Γ ∂a i e 1 + dE Γ, ∂Γ ∂b i e 2 (62) = Γ i -α dE Γ, ∂Γ ∂a i ∂Γ ∂a i (Γ i ) +dE Γ, ∂Γ ∂b i ∂Γ ∂b i (Γ i ) (63)
which is a discrete, constrained evolution [START_REF] Li | Local or global minima: flexible dual front active contours[END_REF]. Noting that ∂Γ/∂a i and ∂Γ/∂b i fit definition (ii) of the pre-defined velocities proposed in Section 3.2.3 and according to the equivalence between discrete and continuous evolutions shown in Section 3.2.4, it can be deduced that the equivalence between ( 26) and (58) holds in the continuous framework. In conclusion, with the polygonal representation, if the pre-defined velocities are chosen equal to the admissible velocities, then the constrained approach is equivalent to the parametric approach.

Uniform cubic B-spline representation If the contour is represented by a uniform cubic B-spline with n segments [START_REF] Brigger | B-spline snakes: A flexible tool for parametric contour detection[END_REF][START_REF] Jacob | A unifying approach and interface for spline-based snakes[END_REF][START_REF] Precioso | Robust real-time segmentation of images and videos using a smooth-spline snake-based algorithm[END_REF][START_REF] Unser | B-spline signal processing: Part I -Theory[END_REF] (see Appendix E for the notations and properties used here), then the parameters involved in the parametric approach are the coordinates (a i , b i ) of each control point q i . Referring to the notations of Section 4.3.1, parameters p are

p = {p j , j ∈ [1, 2n]} = {a 1 , b 1 , a 2 , b 2 . . . , a n , b n } where (a i , b i ) = q i .
(64) For clarity, the following notations will be temporarily used

           dE i a = dE Γ, ∂Γ ∂ai dE i b = dE Γ, ∂Γ ∂bi dE i a,b = ( dE i a dE i b ) T ( 65
)
where A T is the transpose of matrix A. According to (64), procedure (58) can be rewritten as follows

a +1 i = a i -α dE i a b +1 i = b i -α dE i b , i ∈ [1, n] (66)
⇐⇒ q +1 i = q i -α dE i a,b , i ∈ [1, n] . (67)
Spline segment γ i is defined by

γ i (t) = P (t) Q i , t ∈ [0, 1] . (68)
One can verify that, after evolution, the new segment γ +1 i satisfies

γ +1 i (t) = P (t) Q +1 i (69) = γ i (t) -α n j=1 dE Γ, ∂Γ ∂a j ∂γ i ∂a j (t) +dE Γ, ∂Γ ∂b j ∂γ i ∂b j (t) (70)
which is (locally on segment γ i ) the continuous, constrained evolution [START_REF] Kass | Snakes: Active contour models[END_REF]. In conclusion, with the uniform cubic B-spline representation, if the pre-defined velocities are chosen equal to the admissible velocities, then the constrained approach is equivalent to the parametric approach.

About redundancy The admissible velocities are not necessarily mutually orthogonal. In particular, for the polygonal representation, they are not. Yet, these velocities arise in the computation of the gradient of energy (56), a function from R n to R. In this context, the notion of gradient is classical. This might make thinking that redundancy is not always an undesirable property.

Figure 4: Direct/constrained segmentation comparison at a 40-pixel resolution (purposely low to emphasize the difference between the two approaches): initial contour (white plot), result with the direct approach (red plot), with the constrained approach using normal-based pre-defined velocities weighted according to Section 3.2.1 (cyan plot), and with weights computed as suggested in Section 4.2.2 (green plot).

5 Experimental results

Direct vs. constrained approach

A 512×512-pixel image was segmented using the direct then the constrained approach with the same initial contour in a region-competition context (see Fig. 4). Object descriptor (2) was used for both the object and the background and no boundary descriptor was defined. The corresponding energy is

E c (Γ) = Ω (f (x) -µ(Γ)) 2 dx + Ω c (f (x) -µ(Γ c )) 2 dx
(71) and its shape derivative is equal to

dE c (Γ, V ) = - Γ (f (x) -µ(Γ)) 2 -(f (x) -µ(Γ c )) 2 N (s) • V (s) ds . (72)
For the constrained approach, the normal-based predefined velocities (see Section 3.2.3 -definition (i)) were used so that the comparison with the direct approach based by definition on an evolution along the normals (see Section 3.1) is coherent. The active contour was sampled into a polygon with an edge length l of 40 pixels. It was periodically resampled to maintain this resolution.

The algorithm of the direct approach is described in Alg. 1. Constant α min is homogeneous to a number (possibly not an integer) of pixels. It is related to the achievable accuracy of the segmentation (the lower α min , the better the accuracy).

The algorithm of the constrained approach is described in Alg. 2. Note that the weights of the predefined velocities are computed as individual projections of the opposite of the gradient -G, V i . The 40-pixel resolution was purposely low in order to emphasize the difference between the two approaches (which did not allow to segment the object accurately). In fact, in this example, the direct approach segmentation is fairly similar to the segmentation obtained with the constrained approach (see Fig. 4). The constrained approach performed a little bit better, though. Quantitatively, the energy at convergence for the constrained approach (requiring 52 iterations) is 6.5% lower (i.e., better) than the one of the direct approach (requiring 55 iterations). If the weights of the pre-defined velocities are computed so that i β i V i is the projection of -G onto the space of the linear combinations of the pre-defined velocities (see Section 4.2.2), then the energy at convergence (now requiring 71 iterations) is 5.5% lower than the energy reached with the weights based on individual projections. Therefore, when the pre-defined velocities do not form an orthonormal set, the weights should be computed as suggested in Section 4.2.2 for maximum accuracy. In conclusion, the constrained approach can perform better than the direct approach when the contour resolution is low. However, this slight advantage would vanish at high resolution since the error mentioned in Section 3.1.3 would get negligible (see Fig. 5). Experimentally, the segmentation was performed again with a contour resolution of 16 pixels (corresponding to the same ratio with respect to the object size as the one used in Section 5.2). As expected, the results are almost identical for all three methods (see Fig. 6). Therefore, the interest of the constrained approach is not the gain in precision but rather the coherence it provides between the theory developed in the continuous framework and its discrete implementation, and the flexibility it gives to the evolution process, as illustrated in Section 5.2.

Choose an initial polygon Γ = {Γ

i , i ∈ [1, n]} with
an edge length equal to a given resolution l; Choose a threshold α min .

2. Compute domain Ω as the binary mask of Γ.

3. Compute the mean intensities µ(Γ) and µ(Γ c ) in Ω and Ω c , respectively.

4. [Optionally, compute energy E c to check that it decreases.]

5. Compute the L 2 -gradient at Γ i G(Γ, Γ i ) = (f (Γ i )-µ(Γ c )) 2 -(f (Γ i )-µ(Γ)) 2 N (Γ i ) .
(a)

6. Update Γ according to

Γ +1 i = Γ i -α G(Γ, Γ i ) (b)
where α is computed as arg min α≥0 E c (Γ +1 (α)) .

7. If needed, resample Γ to maintain a resolution of l.

8. If α was less than α min , then the algorithm is supposed to have converged; Otherwise go back to step 2.

Algorithm 1: Direct approach algorithm for the polygonal representation.

An example of tracking constraint

An application where the possibility to introduce a priori knowledge in the evolution process can be useful is tracking. Indeed, a usual procedure is to use the segmentation of the object of interest in a frame as the initialization to segment the next frame. Then, the contour deformation necessary to segment the object of interest in the next frame is correlated to the motion of the object6 . Therefore, one can think of choosing motion-based pre-defined velocities. To track a player on the standard test sequence "Football", the normal in the normal-based pre-defined velocities was replaced with a unit vector in the direction of the local motion (see Section 3.2.3 -definition (ii)). The choice of the energy is independent of the choice of the pre-defined velocities. However, it is natural in tracking to use an energy which involves motion [START_REF] Cremers | Variational space-time motion segmentation[END_REF]. Here, considering the complexity of the motion of the object of inter-1.-4. See Alg. 1.

5.

Compute the pre-defined velocities V i according to Section 3.2.3 -definition (i) and compute the shape derivatives dE c (Γ, V i ).

6. Update Γ according to

Γ +1 i = Γ i -α dE c (Γ, V i ) V i (Γ i ) (a) = Γ i -α dE c (Γ, V i ) N (Γ i ) (b)
where α is computed as arg min α≥0 E c (Γ +1 (α)) .

7.-8. See Alg. 1.

Algorithm 2: Constrained approach algorithm for the polygonal representation. est and the slight motion blur, it seemed appropriate to confine the motion use to the pre-defined velocities and to use a motion-free energy able to account for the color variability of the object [2]

E(Γ) = R 2 k(h(Γ, a), h prev (a)) da ( 73
)
where h is a smooth, normalized version of the color histogram in Ω of the frame f t to be segmented

h(Γ, a) = 1 |Ω| Ω g(f t (x) -a) dx ( 74
)
where |Ω| is a measure of Ω and g is a smoothing kernel, e.g., a 2-dimensional Gaussian. Similarly, h prev is a smooth, normalized version of the color histogram of the segmentation in the previously segmented frame f t-1 . Note that a and f t (x) should belong to R 3 . However, to limit the computation load, only the two most significant color components out of Figure 6: Direct/constrained segmentation comparison at a 16-pixel resolution: result with the direct approach (red plot), with the constrained approach using normal-based pre-defined velocities weighted according to Section 3.2.1 (cyan plot), and with weights computed as suggested in Section 4.2.2 (green plot). three were considered. Function k is a positive function from R 2 to R defined as k(x, y) = (x -y) 2 . A maximal area constraint [START_REF] Roy | Segmentation of a vector field: Dominant parameter and shape optimization[END_REF] was added to energy (73) since its sensitivity decreases in the inner neighborhood of the correct segmentation

E A (Γ) = -δ Ω dx ( 75
)
where δ is a weighting parameter. The shape derivative of the sum of ( 73) and ( 75) is

dE(Γ, V ) = - Γ 1 |Ω| (g ⋆ k ′ (h(Γ, •), h ref (•)))(f t (s)) - 1 |Ω| 2 R 2 h(Γ, a) k ′ (h(Γ, a), h ref (a)) da -δ N (s) • V (s) ds ( 76
)
where ⋆ is the convolution operator and k ′ is the partial derivative of k with respect to its first variable. The minimization algorithm follows the same steps as Alg. 2. Resolution l was chosen equal to 10 pixels (to be compared to a size of frame of 352×288 pixels). In step 5 of Alg. 2, the local motion at the vertices, required to determine the pre-defined velocities, was estimated by a sub-optimal block matching technique [START_REF] Zhu | A new diamond search algorithm for fast block matching motion estimation[END_REF] with a 1/4-pixel precision. Here, a block was defined as the intersection between a 21×21-pixel square centered on a vertex and the mask of the current active contour. This prevented the pixels considered to be outside the object from interfering with the motion estimation. This procedure is given for illustrative purposes only.

In a real-world tracking application, motion estimation should certainly be more sophisticated [START_REF] Arulampalam | A tutorial on particle filters for online nonlinear/non-Gaussian Bayesian tracking[END_REF][START_REF] Isard | CONDENSATION -Conditional density propagation for visual tracking[END_REF][START_REF] Robinson | Fast local and global projection-based methods for affine motion estimation[END_REF]. Figure 7 presents the results obtained with the constrained approach using normal-based pre-defined velocities and motion-based pre-defined velocities. The normal-based segmentation is globally similar to the motion-based segmentation except for a small region wrongly included above the player. These segmentations correspond to two local minima, one of which being more relevant. In a way, the motion-based evolution took the shortest path toward the object of interest and converged toward a more satisfying minimizer. As a matter of fact, the convergence was reached after 10 iterations while it took 14 iterations for the normal-based evolution. An intuition of the behavior of the normal-based evolution is given by the normals shown in Fig. 7 compared to the local motions. When a normal has a direction close to the direction or opposite direction of the local motion7 , then both evolution processes behave similarly. However, when these directions are close to be orthogonal, then the normal-based process could have a tendency to evolve toward a less relevant local minimum and/or to require more iterations to converge. Indeed, the local motion is likely to be a better direction to take in a tracking application.

An example of symmetry constraint

Pre-defined velocities can also be chosen in order to maintain symmetry. Suppose that the initial contour has an axis of symmetry. Although the previous examples used one pre-defined velocity per sample, it is not necessary to do so. A pre-defined velocity could be assigned to each pair of symmetric samples, the velocity being itself symmetric with respect to the axis. For example, if the contour is represented by a polygon, then the sum for a pair of symmetric samples of the normalbased pre-defined velocities proposed in Section 3.2.3definition (i) could specify one such pre-defined velocity. The weight of this velocity being unique, the displacements of the two related samples will be identical in norm. Since their directions is given by the local normal, which is symmetric with respect to the axis of symmetry, the symmetry of the contour will be maintained during evolution.

To illustrate this procedure, an object artificially made slightly asymmetric was segmented using energy (71) and Alg. 2 (see Fig. 8). In this simple example, the axis of symmetry was supposed to be known. (For comparison purposes, the segmentation was also computed using one normal-based pre-defined velocity per sample, as in Section 5.1.) Note that this method does not allow to recover the symmetric object. The computed segmentation is rather a symmetric shape representing a trade-off between asymmetric regions kept inside and left outside the contour.

Conclusion

A variational approach to segmentation was presented and the shape derivative was highlighted as a convenient framework for deriving the evolution equation of an active contour from the energy to be minimized. It was noted that an error is induced by a discretization Figure 8: Segmentation of an asymmetric object with (green contour) and without (red contour) a symmetry constraint imposed by an appropriate choice of predefined velocities.

of the evolution equation when the contour velocity is taken equal to the gradient of the energy associated with the L 2 inner product. Although there are no consequences in practice if the contour is discretized finely enough, the constrained approach proposed to solve this problem also gives more flexibility to the active contour process by allowing to introduce some a priori knowledge. This property can be useful, e.g., for a tracking application where motion gives an a priori on the position of the tracked object in the next frame.

A Shape gradient of a boundary energy

Let us consider the following boundary energy

E(Γ) = Γ ϕ(s) ds ( 77
)
where Γ is the oriented boundary ∂Ω of an open set Ω of R 2 and s is the arc-length parameterization of Γ. The shape derivative of (77) is equal to [START_REF] Sokolowski | Introduction to shape optimization: Shape sensitivity analysis[END_REF][START_REF] Delfour | Shapes and geometries: Analysis, differential calculus and optimization[END_REF]

dE(Γ, V ) = Γ ∂ϕ(s) ∂N -ϕ(s) κ(s) N (s) • V (s) ds
(78) where N is the inward unit normal of Γ and κ is the curvature of Γ.

Note that this result can also be obtained by a calculus of variations [START_REF] Caselles | Geodesic active contours[END_REF].

B Shape gradient of a domain energy

Let Ω be an open set of R 2 and let T be a transformation of Ω such that

   Ω = T (τ = 0, Ω) Ω T (τ ) = T (τ, Ω) x T (τ ) = T (τ, x), x ∈ Ω . ( 79
)
The deformation at x is defined as

V T (x) := lim τ →0 x T (τ ) -x τ (80) = ∂T ∂τ (0, x) . (81)
For clarity, Ω T (τ ), x T (τ ), and V T (x) are referred to as Ω(τ ), x(τ ), and V (x), respectively. Function V is called the velocity of Ω.

Let us consider the following domain energy

E(Γ) = Ω φ(Γ, x) dx ( 82
)
where Γ is the oriented boundary ∂Ω of Ω. The same energy computed for the transformed domain Ω(τ ) is equal to

E(Γ, T, τ ) = Ω(τ ) φ(Γ(τ ), x) dx . (83)
Then, the shape derivative of (82) is defined as

dE(Γ, T ) := lim τ →0 E(Γ, T, τ ) -E(Γ, T, 0) τ . (84)
It is equal to [START_REF] Sokolowski | Introduction to shape optimization: Shape sensitivity analysis[END_REF][START_REF] Delfour | Shapes and geometries: Analysis, differential calculus and optimization[END_REF]

dE(Γ, T ) = Ω ∂φ(Γ(τ ), x) ∂τ τ =0 dx - Γ φ(Γ, s) N (s) • V (s) ds ( 85
)
where s is the arc-length parameterization of Γ and N is the inward unit normal of Γ. Since V appears explicitly in the shape derivative expression, dE(Γ, T ) can be replaced with dE(Γ, V ). Note that this result can also be obtained by a calculus of variations [START_REF] Aubert | Image segmentation using active contours: Calculus of variations or shape gradients?[END_REF].

C Rewriting the shape derivative as a boundary integral

Under weak assumptions, the shape derivative of a domain energy (see Appendix B) has an equivalent expression in the form of a boundary integral [START_REF] Sokolowski | Introduction to shape optimization: Shape sensitivity analysis[END_REF][START_REF] Solem | A geometric formulation of gradient descent for variational problems with moving surfaces[END_REF]. Such an expression is convenient in the active contour framework (see Section 3.1.1). Indeed, the evolution equation of an active contour can be easily deduced from a boundary integral-only expression. Two (sets of) conditions are proposed to allow to practically rewrite shape derivative (85) into an expression without any domain integral.

C.1 Recursive applications of the shape derivative

Under some conditions on the dependency of φ on Γ, applying recursively (85) to its first integral leads to an expression containing no domain integral [START_REF] Jehan-Besson | DREAM 2 S: Deformable regions driven by an eulerian accurate minimization method for image and video segmentation[END_REF] dE(Γ,

V ) = - Γ Ψ(Γ, s) N (s) • V (s) ds . ( 86
)
Let us assume that

   φ i (Γ, x) = φ i (g i (Γ), x), i ∈ [1, n -1] g i (Γ) = Ω φ i+1 (Γ, x) dx, i ∈ [1, n -1] φ n (Γ, x) = φ n (x) (87)
where φ 1 = φ (note that g i has the same form as ( 82)).

The first integral of (85) reads

Ω ∂φ 1 (g 1 (Γ(τ )), x) ∂τ τ =0 dx = Ω ∂φ 1 (k, x) ∂k g1(Γ) ∂g 1 (Γ(τ )) ∂τ τ =0 dx (88) = ∂g 1 (Γ(τ )) ∂τ τ =0 Ω ∂φ 1 (k, x) ∂k g1(Γ) dx (89) := dg 1 (Γ, V ) Ω ∂φ 1 (k, x) ∂k g1(Γ) dx (90) = dg 1 (Γ, V ) A 1 (Γ) (91)
where dg 1 (Γ, V ) is equal to

dg 1 (Γ, V ) = Ω ∂φ 2 (Γ(τ ), x) ∂τ τ =0 dx - Γ φ 2 (Γ, s) N (s) • V (s) ds . (92)
Then, the development leading to (91) can be repeated with the successive domain integrals present in dg i (Γ, V ), i increasing, until

dg n-1 (Γ, V ) = - Γ φ n (s) N (s) • V (s) ds (93)
which does not contain any domain integral since φ n is independent of Γ. Gathering the successive boundary integrals together, the shape derivative is equal to

dE(Γ, V ) = - Γ Ψ(Γ, s) N (s) • V (s) ds (94) = - Γ   n i=1 φ i (g i (Γ), s) i-1 j=1 A j (Γ)   N (s) • V (s) ds ( 95
)
where A j (Γ) is equal to

A j (Γ) = Ω ∂φ j (x, k) ∂k gj (Γ) dx ( 96
)
with the following convention

φ n (g n (Γ), s) = φ n (s) . (97)
For instance, φ may involve the variance g 1 of f in Ω. The variance involves the average value g 2 of f in Ω. Finally, the average value involves f but no term depending on Γ.

C.2 Domain integral equal to zero

Under some conditions, Ψ in (86) is simply equal to φ [START_REF] Roy | Segmentation of a vector field: Dominant parameter and shape optimization[END_REF]. In other words, the first integral of (85) is equal to zero. For instance, if φ is given by (2), i.e., (f (x)-µ(Γ)) 2 , the integrand of the first integral of (85) is equal to

∂φ(Γ(τ ), x) ∂τ τ =0 = ∂ (f (x) -µ(Γ(τ ))) 2 ∂τ τ =0 (98) = -2 (f (x) -µ(Γ)) ∂µ(Γ(τ )) ∂τ τ =0 (99)
where µ(Γ) can be written as follows

µ(Γ) = Ω f (x) dx / Ω dx . (100)
Then, the first integral of (85) is equal to

Ω ∂φ(Γ(τ ), x) ∂τ τ =0 dx = -2 ∂µ(Γ(τ )) ∂τ τ =0 Ω (f (x) -µ(Γ)) dx (101) = -2 ∂µ(Γ(τ )) ∂τ τ =0 Ω f (x) dx -µ(Γ) Ω dx (102) = 0 . (103)
In general, a sufficient condition for the first integral of (85) to be equal to zero is

φ(Γ, x) = φ(g(Γ), x) g(Γ) = arg min k Ω φ(k, x) dx . (104)
In other words, g(Γ) is the minimizer of (82) seen as a function of g with Γ fixed. For convenience, the following notation will be used

g(Γ) = arg min k E Γ (k) . ( 105
)
If g is assumed to be differentiable, the development of (91) can be continued using the new condition ( 104)

Ω ∂φ(g(Γ(τ )), x) ∂τ τ =0 dx = dg(Γ, V ) Ω ∂φ(k, x) ∂k g(Γ) dx (106) = dg(Γ, V ) ∂ ∂k Ω φ(k, x) dx g(Γ) (107)
= dg(Γ, V ) ∂E Γ (k) ∂k g(Γ)
.

If g(Γ) is a constraint-free minimizer of E Γ , then it can be immediately concluded that (108) is equal to zero. Otherwise, let us assume that g(Γ) is the minimizer of E Γ under the constraint ξ(g(Γ)) = 0 .

Then, there exists a Lagrange multiplier λ such that

∂E Γ (k) ∂k g(Γ) = λ ∂ξ(k) ∂k g(Γ) . (110)
Therefore,

Ω ∂φ(g(Γ(τ )), x) ∂τ τ =0 dx = λ ∂g(Γ(τ )) ∂τ τ =0 ∂ξ(k) ∂k g(Γ(τ =0)) (111) = λ ∂ξ(g(Γ(τ ))) ∂τ τ =0 . ( 112
)
By definition, for any τ , ξ(g(Γ(τ ))) has the same value (equal to zero). It can be concluded that (112) is equal to zero.

As suggested by a reviewer, noting that E(Γ) can be written as F (Γ, g(Γ)) leads to an immediate proof of the result of C.2.

D Validity of the direct approach

As suggested by a reviewer, despite the discretization error, the direct approach can still be valid (i.e., evolution ( 16) generates a minimizing sequence) if the contour discretization is fine enough. The following development is largely inspired from the reviewer's proof.

Let G be the gradient associated with the

L 2 inner product G(Γ, s) = -Ψ(Γ, s) N (s) . (113)
For convenience, G(s) and Ψ(s) will be used instead of G(Γ, s) and Ψ(Γ, s), respectively. Let us assume that Γ is represented by a polygon with n edges of equal length l (see Section 1.4 for notations). Note that l is equal to s i+1 -s i . Let L be the length of Γ. Let N be a C 2 pseudo-normal of Γ. Function Ψ is also assumed to be C 2 so that G ′′ exists and is continuous. With the direct approach, the induced velocity Ṽ at s ∈

[s i , s i+1 ] is equal to Ṽ (s) = -(1 -s) G(s i ) -s G(s i+1 ) (114) = Ψ(s i ) (1 -s) N (s i ) +Ψ(s i+1 ) s N (s i+1 ) (115)
where s is equal to (s -s i )/l. Note that (1 -s) N (s i ) and s N (s i+1 ) are equal to V i (s) and V i+1 (s), respectively, the normal-based pre-defined velocities proposed in Section 3.2.3 -definition (i). Also note that V j (s) is equal to zero if j is not equal to i or i + 1. Therefore, the induced velocity can be rewritten as follows

Ṽ (s) = n i=1 Ψ(s i ) V i (s) (116)
where s can now belong to [0, L]. The shape derivative (9) at Ṽ is equal to

dE(Γ, Ṽ ) = Γ G(s) • n i=1 Ψ(s i ) V i (s) ds (117) = n i=1 Ψ(s i ) Γ G(s) • V i (s) ds . ( 118
)
Velocity V i being different from zero on γ i-1 and γ i only, the shape derivative is equal to

dE(Γ, Ṽ ) = n i=1 Ψ(s i ) γi-1 G(s) • V i (s) ds + γi G(s) • V i (s) ds . (119)
According to Taylor's theorem with the Lagrange form of the remainder term, there exists

t i in ]s i-1 , s i+1 [ such that G(s) = G(s i )+G ′ (s i ) (s-s i )+ G ′′ (t i ) 2 (s-s i ) 2 . ( 120
)
The first integral of (119) is equal to

I 1 = G(s i ) • γi-1 V i (s) ds + G ′ (s i ) • γi-1 (s -s i ) V i (s) ds + G ′′ (t i ) 2 • γi-1 (s -s i ) 2 V i (s) ds . (121)
With the change of variable t = 2s i -s, noting that 2s i -s i-1 is equal to s i+1 , and due to the symmetry of V i with respect to s i , one can verify that γi-1

(s -s i ) V i (s) ds = - γi (s -s i ) V i (s) ds . (122) Therefore, (119) is equal to dE(Γ, Ṽ ) = n i=1 Ψ(s i ) G(s i ) • Γ V i (s) ds + G ′′ (t i ) 2 • Γ (s -s i ) 2 V i (s) ds . (123)
Replacing V i with its expression, one can verify that Γ V i (s)ds is equal to l N (s i ) and that

Γ (s -s i ) 2 V i (s) ds = l 2 6 Γ V i (s) ds . (124)
Finally, (119) is equal to

dE(Γ, Ṽ ) = n i=1 Ψ(s i ) G(s i ) + l 2 12 G ′′ (t i ) • Γ V i (s) ds (125) = n i=1 Ψ(s i ) (-Ψ(s i ) N (s i ) + l 2 12 G ′′ (t i ) • l N (s i ) (126) = l n i=1 Ψ(s i ) l 2 12 G ′′ (t i ) • N (s i ) -Ψ(s i ) 2 . ( 127
)
If l is such that

l 2 n i=1 Ψ(s i ) 12 G ′′ (t i ) • N (s i ) ≤ n i=1 Ψ(s i ) 2 , (128)
then the negativity of (127) is guaranteed. The following stronger condition can be deduced from (128)

l 2 M 12 n i=1 |Ψ(s i )| ≤ n i=1 Ψ(s i ) 2 (129) ⇐⇒ l 2 ≤ 12 M n i=1 Ψ(s i ) 2 n i=1 |Ψ(s i )| (130) ⇐⇒ l 2 ≤ 12 M l n i=1 Ψ(s i ) 2 l n i=1 |Ψ(s i )| (131)
where M is equal to

M = sup s |G ′′ (s)| (132)
and M and n i=1 |Ψ(s i )| are assumed to be different from zero (otherwise, the negativity of (127) is trivially verified). Functions |Ψ| and Ψ 2 being Riemannintegrable, the limit when l tends toward zero of the right-hand side of (131) exists and is equal to 812 M Γ Ψ(s) 2 ds Γ |Ψ(s)| ds

, ( 133
)
assuming that Ψ is not identically equal to zero on Γ (otherwise, the negativity of ( 127) is again trivially verified). Limit (133) being strictly positive, there exists a value l max such that if l is smaller than l max , then condition (131) is verified, implying the negativity of (127).

Let I u and S u be the integral and the sum at the numerators of ( 133) and (131), respectively. The error between I u and S u is such that where l root max is the largest root of the polynom in (139). One can check that it exists and is positive.

Note that l root max depends on Γ and Ψ, which implies that it should be computed at each iteration of the active contour evolution.

The reviewer's proof led to the following, much simpler upper bound 1/ M √ L.

E Some properties of closed, uniform cubic B-splines

Here, we present the properties of closed, uniform cubic B-splines that are used in Sections 3.2.4 and 4.3.2.

A closed, uniform cubic B-spline with n segments [START_REF] Brigger | B-spline snakes: A flexible tool for parametric contour detection[END_REF][START_REF] Jacob | A unifying approach and interface for spline-based snakes[END_REF][START_REF] Precioso | Robust real-time segmentation of images and videos using a smooth-spline snake-based algorithm[END_REF][START_REF] Unser | B-spline signal processing: Part I -Theory[END_REF] is uniquely defined by the coordinates (a i , b i ) of n control points q i and the following blending function B (see Fig. 9)

B(t) =              2 3 -t 2 + |t| 3 2 if 0 ≤ |t| < 1 (2 -|t|) 3 6 if 1 ≤ |t| < 2 0 if |t| ≥ 2 . ( 141
)
The equation of spline segment γ i is , i.e., blending function B shifted to be "centered" on Γ i and "seen" from segment γ i+1 .

γ i (t) = i+2 j=i-1 q j B i j (t), t ∈ [0, 1] (142)
where B i j represents blending function B shifted to be "centered" on Γ j and "seen" from segment γ i (see Fig. 9). Mathematically,

B i j (t) = B(t -(j -i)) . (143)
Note that the following notation abuses are made    q 0 = q n q n+1 = q 1 q n+2 = q 2 .

(144)

The extremities γ i (0) and γ i (1) of a segment are equal to sampling points Γ i and Γ i+1 , respectively.

Property 1 Although B i j (t) was defined for t ∈ [0, 1], the domain of definition can be extended to [-n/2, n/2], covering the entire curve. Note that if t is an integer in this extended interval, then it corresponds to Γ i+t .

Property 2 Sampling points and control points are linked together by the following equation

Γ i = 1 6 (q i-1 + 4q i + q i+1 ) . ( 145
)
Property 3 Replacing B with its definition, the segment equation reads γ i (t) = t 3 6 (q i+2 -3q i+1 + 3q i -q i-1 ) + t 2 2 (q i+1 -2q i + q i-1 ) + t 2 (q i+1 -q i-1 ) + 1 6 (q i+1 + 4q i + q i-1 ), Q i = ( q i+2 q i+1 q i q i-1 ) T (148) and A T is the transpose of matrix A. Let (e 1 , e 2 ) be the canonical basis of R 2 . Note that for any set of α j 's, and similarly when taking the derivative with respect to b j .

P (t) = ∂γ i ∂a i+2 (

Figure 1 :

 1 Figure 1: Incorrect deformation due to the sampling of V . (a) Disks: polygon vertices; Solid line: polygon before deformation; Dashed line: polygon defined by the translated vertices; Dotted line: polygon deformed according to V . (b) Disks: curve samples; Solid line: curve before deformation; Dashed line: curve interpolating the translated samples; Dotted line: curve deformed according to V .

Figure 2 :

 2 Figure 2: A possible choice for the pre-defined velocity V i .

Figure 3 :

 3 Figure 3: Transformation of point x on edge γ i . tion 3.2.3 (either of the two definitions). The velocity V (t) at γ i (t) is the combination of two velocities

Figure 5 :

 5 Figure 5: The error between the correctly deformed curve (dotted line) and the wrongly estimated curve (dashed line) decreases when resolution gets higher. Disks: curve samples; Solid line: curve before deformation; Dashed line: curve interpolating the (correctly) translated samples; Dotted line: curve correctly deformed.

Figure 7 :

 7 Figure 7: Tracking of a player on frame 75 of the standard test sequence "Football" with the constrained approach. Top left: manual segmentation of frame 74 (used as the initial contour for the segmentation of frame 75); Top right: histogram (for display purposes, upper and lower zero-valued regions were cut out) of the manually segmented region (used as histogram h prev ); Middle left: initial contour on frame 75 (copied from frame 74); Middle right: a closeup of the normals (red arrows) and the local motions (green arrows) computed at the vertices of the polygon at the first iteration; Bottom left: segmentation of frame 75 using normal-based pre-defined velocities; Bottom right: segmentation of frame 75 using motion-based pre-defined velocities.

  |I u -S u | ≤ L P l (134) =⇒ I u -L P l ≤ S u (135) where P is equal to 2 sup s |Ψ(s) Ψ ′ (s)|. Let I d and S d be the integral and the sum at the denominators of (133) and (131), respectively. Using the Chasles's equality and the mean value theorem, one can check that S d -I d = l n i=1 |Ψ(s i )| -|Ψ(t i )| (136) ≤ n l (sup s |Ψ(s)| -inf s |Ψ(s)|) Q (137) where t i belongs to ]s i , s i+1 [. Then, a condition stronger than (131) is l 2 ≤ 12 M I u -L P l I d + L Q (138) ⇐⇒ M 12 (I d + L Q) l 2 + L P l -I u ≤ 0 (139) ⇐⇒ l ≤ l root max (140)

Figure 9 :

 9 Figure 9: Top: the uniform cubic B-spline blending function B. It is equal to zero outside the interval [-2, 2]; Bottom: function B i+1 i

  where α is a polynom. If j is outside the interval [i -1, i + 2], then ∂γ i

										∂a j	(t) = 0	(152)
									since, in this case, γ i is independent of a j . In both
									cases,	∂γ i ∂a j	(t) • e 2 = 0 .	(153)
									Similarly,	∂γ i ∂b j	(t) • e 1 = 0 .	(154)
									Property 5 Considering (151) and (152), it can be
									deduced that
									i+2 j=i-1	α j	∂γ i ∂a j	(t) =	n j=1	α j	∂γ i ∂a j	(t)	(155)
				t) • e 1	∂γ i ∂a i+1	(t) • e 1
		∂γ i ∂a i	(t) • e 1	∂γ i ∂a i-1	(t) • e 1	(149)
	=	∂γ i ∂b i+2	(t) • e 2	∂γ i ∂b i+1	(t) • e 2
		∂γ i ∂b i	(t) • e 2	∂γ i ∂b i-1	(t) • e 2 . (150)

Property 4 Let j be an integer in the interval [i -1, i + 2]. Then, the derivative of (146) with respect to a j is equal to

∂γ i ∂a j (t) = α(t) e 1

(151)

For convenience, φ and ϕ will be used instead of φ f and ϕ f .

A problem of maximization is trivially turned into a problem of minimization.

The sequence x(n + 1) = f (x(n)), x(0) = x 0 is denoted by x +1 = f (x), x 0 .

In particular,[START_REF] Debreuve | Active contour segmentation with a parametric shape prior: link with the shape gradient[END_REF] does certainly not transform, in general, a polygon into another polygon.

For definitions and notations, see Appendix B.

Ideally, the motion could even transform the segmentation in a frame directly into the segmentation in the next frame. However, because the sequence is a two-dimensional projection of a three-dimensional scene and because the available motion is usually an apparent motion, this is not the case in practice.

If the direction of the normal is opposite to the direction of the local motion, the shape gradient takes equal values in absolute value for both but with opposite signs, making no difference as far as evolution is concerned.

Remember that sn is equal to L -l (see Section 1.4). Therefore, it is correct to sum for all i up to n (and not up to n -1) to approximate the integral over Γ.

Acknowledgments

The authors would like to thank Jean-Paul Marmorat, École des Mines de Paris, C.M.A. Sophia Antipolis, France, for a fruitful discussion on the shape derivative in the context of a parametric active contour.