Eric Debreuve
email: ede-breuv@doug.med.utah.edu

Member, IEEE Michel Barlaud

Gilles Aubert

Member, IEEE Ivan Laurette

Member, IEEE Jacques Darcourt

Éric Debreuve

§ Preprint

§ G Aubert

Space-time segmentation using level set active contours applied to myocardial gated SPECT

Keywords: 3-D gated SPECT sequence, space-time segmentation, geometrical constraint, active contour

Introduction

Classically, when segmenting a sequence, each frame of the sequence is segmented independently from the others [START_REF] Germano | Automatic quantification of ejection fraction from gated myocardial perfusion SPECT[END_REF]- [START_REF] Faber | Left ventricular function and perfusion from gated SPECT perfusion images: An integrated method[END_REF]. Alternatively, each frame can be segmented taking into account other previous and subsequent frames, or even the entire sequence. This second approach, called space-time segmentation, offers the advantage of providing a time consistent segmentation.

Kornprobst et al. [START_REF] Kornprobst | Image sequence analysis via partial differential equations[END_REF] developed a variational method for space-time segmentation in the scope of video monitoring. The method assumes that such sequences are composed of moving (and warping) objects against a still background (in other words, the frames are all registered with respect to one another). The aim is to detect the moving objects in the sequence, i.e. to segment the objects in each frame of the sequence. Their criterion takes into account the entire sequence causing individual frame segmentations to depend on all the others. Nevertheless, it is geometrically unconstrained. As a consequence, the segmentation is performed pixel by pixel.

Pixel based segmentation methods should be distinguished from region based segmentation methods. In the case of sequence segmentation, a region based method is developed in [START_REF] Paragios | Geodesic active regions for motion estimation and tracking[END_REF]. Assuming that a background reference image is known, the segmentation is performed using active contours whose evolution velocity depends on intensity and motion information. Part of the motion information is modeled utilizing the optical flow constraint. This approach is statistical in the sense that intensity values are widely replaced with probability density functions. The authors discuss the reliability of the different evolution terms. Another region based statistical approach is presented in [START_REF] Yezzi | A statistical approach to snakes for bimodal and trimodal imagery[END_REF] for single image segmentation. It also uses active contours. The evolution velocity is such that it increases the discrepancy between statistical values computed in each region, or equivalently, such that it maximizes the homogeneity in each region. The number of independent statistical values required is directly related to the number of types of region, which must be chosen ini-tially. The active contour process is regularized using an additional curvature term.

Based on the work in [START_REF] Kornprobst | Image sequence analysis via partial differential equations[END_REF], we propose to geometrically constrain the segmentation process to make it region based. In other words, instead of looking for a given pixel/voxel to belong to an object or to the background, our goal is to search the set of pixels/voxels that form the same object; The remaining pixels/voxels belong to the background. To achieve this goal, we have developed a new space-time segmentation criterion (restricted to a unique moving and warping object) which involves the space domains of the object in each frame. The derivation of the Euler-Lagrange equations is constrained by the object domain boundaries. As a consequence, we find an expression of an evolution velocity normal to the object boundary. This naturally implies an active contour process [START_REF] Kass | Snakes: Active contour models[END_REF]- [START_REF] Caselles | Geodesic active contours[END_REF] (further theoretical studies and applications can be found in [START_REF] Yezzi | A geometric snake model for segmentation of medical imagery[END_REF], [START_REF] Siddiqi | Area and length minimizing flows for shape segmentation[END_REF]- [START_REF] Angenent | On the Laplace-Beltrami operator and brain surface flattening[END_REF]). We have implemented it using a level set modelization [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF]- [START_REF] Sethian | Level set methods[END_REF] which guarantees that the segmentation contours (the zero levels) are closed hyper-surfaces [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF], [START_REF] Debreuve | Attenuation map segmentation without reconstruction using a level set method in nuclear medicine imaging[END_REF].

We have tested this space-time segmentation method on simulated and clinical 3-D myocardial gated SPECT sequences, for which dynamic features such as myocardial wall thickening or ejection fraction may be computed from the endocardial and epicardial limit segmentation in each frame of the sequence. As an illustration, we have computed the ejection fractions of the segmented sequences using the sign of level sets.

Since the method developed in [START_REF] Kornprobst | Image sequence analysis via partial differential equations[END_REF] is the starting point of our work, we briefly sketch some of its key components in section 3. Section 4 presents the development of our space-time segmentation algorithm. In section 5 we present results on simulated and clinical data. Finally, we discuss relevant issues in section 6.

Notations

Two-dimensional or three-dimensional images are denoted with an additional "(.)" or "(., .)". t is the iteration index when embedding an unknown into a family of evolving terms converging toward the solution. For example, unknown a or a(.) may be denoted a(t) or a(., t) respectively when considering the iterative resolution process. Vectors are denoted with an overhead arrow. Tables 1 and2 present the main symbols used in this paper. 3 Pixel based segmentation method

Unconstrained criterion

Let us assume that a sequence is composed of frames with a still background occluded by some moving (and warping) objects (for example, a sequence acquired by a static video monitoring camera). Therefore, the frames are all registered with respect to one another so that a particular background pixel, when not occluded by an object, can be assumed to be of fixed intensity. Then, if the background image is known (a frame where the moving objects have been removed, uncovering the background), it can be compared pixel by pixel to any frame of the sequence. If two corresponding pixel values are equal, the pixel belongs to the background. If two corresponding pixel values are different, the pixel belongs to a moving object. Based on these assumptions and properties, the following criterion was developed for space-time segmentation in the scope of video monitoring [START_REF] Kornprobst | Image sequence analysis via partial differential equations[END_REF]:

J 1 (c n , b) = N n=1 Ω c 2 n (x) × {b(x) -S n (x)} 2 +α Ω {c n (x) -1} 2 +β Ω ϕ b (|∇b(x)|) + γ N n=1 Ω ϕ c (|∇c n (x)|) (1)
where S n (.) is frame n of the sequence to segment, c n (.) is the unknown space-time segmentation (see the

(x) = N n=1 ĉ2 n (x)S n (x) N n=1 ĉ2 n (x) , x ∈ Ω (2)
where ĉn (x) is such that:

• ĉn is a function from Ω to [0, 1],

• ĉn (x) is close to one if x belongs to the background in frame n, and

• ĉn (x) is close to zero if x belongs to a moving object in frame n.

Without the ϕ c -segmentation regularization term (γ = 0), ĉn (x) equals:

ĉn (x) = α α + {b(x) -S n (x)} 2 • (3)
Let us give an intuitive explanation of criterion [START_REF] Germano | Automatic quantification of ejection fraction from gated myocardial perfusion SPECT[END_REF]. The first sum over n is the data consistency term. The last two terms are dedicated to the regularization of the solutions. The data consistency term is composed of two integrals of positive expressions. In order to minimize criterion (1), these expressions must be equal to zero or close to zero. The first integral equals zero if and only if for all x: b(x) = S n (x) or c n (x) = 0 or both. If b(x) = S n (x) and c n (x) = 0, then the expression in the second integral is maximum (since: c n (x) ∈ [0, 1]). Therefore, the two conditions b(x) = S n (x) and c n (x) = 0 seem mutually exclusive. Moreover, when b(x) = S n (x), the expression of the first integral equals zero regardless of what the value of c n (x) may be. This value should be one so that the expression of the second integral equals zero. In summary, there are two cases that allow minimization of the criterion:

b(x) = S n (x) c n (x) = 1 or b(x) = S n (x) c n (x) = 0 • (4)
We can then conclude that a point x of the background in frame n is characterized by c n (x) = 1 (because the value S n (x) of the frame is equal to the background value b(x)). Conversely, a point x of a moving object in frame n is characterized by c n (x) = 0. In this second case, the expression in the second integral has the maximum value. As a consequence, c n (x) converges toward a small value in order to minimize the first integral without maximizing the second one. This value of convergence depends on parameter α (it is an increasing function of α). Therefore, the second integral acts as a penalty for assigning a particular point to a moving object. Without this penalty criterion (1) would be minimized with c n (x) = 0 for all n and all x, i.e. each frame is entirely composed of moving objects.

This segmentation model has a long history in the literature. In its earliest form (in the context of segmenting single images), it goes back to the variational formulation of Mumford and Shah [START_REF] Mumford | Boundary detection by minimizing functionals[END_REF]. Let us try to see why. We consider the single image case (N = 1). Let us assume that c 1 (.) is the characteristic function over b(.), i.e. c 1 (x) equals one if x belongs to the background and otherwise equals zero. The Mumford-Shah criterion for image segmentation is:

J M-S (w, Γ) = Ω\Γ {w(x) -S 1 (x)} 2 +γ Ω\Γ |∇w(x)| 2 + L (5)
where w(.) is the segmented image, Γ represents the contours in S 1 (.), and L is a measure of the contour lengths. The pair (w(.), Γ) represents the unknown segmentation of S 1 (.). We still consider the case of objects on a background. Therefore, w(x) can be written:

c 1 (x)b(x) + (1 -c 1 (x))o 1 (x)
where o 1 (x) denotes the object image. Then, the term in the first integral of (5) (Let us denote: A(x) = {w(x) -S 1 (x)} 2 ) equals:

A(x) = {[c 1 (x)b(x) + (1 -c 1 (x))o 1 (x)] -[c 1 (x)S 1 (x) + (1 -c 1 (x))S 1 (x)]} 2 (6) = c 1 (x) 2 {b(x) -S 1 (x)} 2 + (1 -c 1 (x)) 2 {o 1 (x) -S 1 (x)} 2 + 2c 1 (x)(1 -c 1 (x)){b(x) -S 1 (x)} {o 1 (x) -S 1 (x)} • (7)
Expression [START_REF] Faber | Left ventricular function and perfusion from gated SPECT perfusion images: An integrated method[END_REF] has similarities with criterion [START_REF] Germano | Automatic quantification of ejection fraction from gated myocardial perfusion SPECT[END_REF]. In particular, the first term of (7) can be found in the first integral of (1). The second term of (7) represents the term in the second integral of (1) (the exact quadratic error {o 1 (x) -S 1 (x)} 2 has been replaced with parameter α that serves as a mean quadratic error). The third term of (7) equals zero since c 1 (x) equals either zero or one. The second term of (5) corresponds to the ϕ c -regularization term of (1).

As a conclusion about the method of Kornprobst et al. [START_REF] Kornprobst | Image sequence analysis via partial differential equations[END_REF], each frame can be segmented into background and objects using a threshold value between zero and one. Every pixel with a c n value below this threshold corresponds to an object. Every pixel with a c n value above this threshold corresponds to the background.

Using this criterion, the computed segmentation is not globally constrained: There are no assumptions that compel some pixels to belong to the same moving object. Regardless of the ϕ c -regularization term, the segmentation process is performed pixel by pixel, and is therefore unconstrained.

Space-time segmentation algorithm

The space-time segmentation is represented by the sequence ĉn (.) which results from the minimization of criterion [START_REF] Germano | Automatic quantification of ejection fraction from gated myocardial perfusion SPECT[END_REF]. The minimization algorithm developed in [START_REF] Kornprobst | Image sequence analysis via partial differential equations[END_REF] is based on alternate resolutions of the Euler-Lagrange equations with respect to b(.) and c n (.). For both b(.) and c n (.), these equations are partial differential equations involving a weighted Laplacian [START_REF] Geman | Nonlinear image recovery with half-quadratic regularization[END_REF], [START_REF] Charbonnier | Deterministic edge-preserving regularization in computed imaging[END_REF]. They are discretized and solved using a Gauss-Seidel iterative algorithm.

4 New region based segmentation method 4.1 Object/background balanced criterion

Criterion (1) gives priority to the background with respect to the moving objects. Indeed, the term Ω c 2 n (x) × {b(x) -S n (x)} 2 involves the background image b(.) whereas the term α Ω {c n (x) -1} 2 does not explicitly involve the values of the moving objects. The reason is that moving objects in video monitoring sequences are generally not homogeneous. On the contrary, in SPECT images, the observed organ pixels/voxels are roughly of the same value. Then, a sequence of this type of images can be considered as one or several homogeneous moving (and warping) organs on a still background with all of the images registered with respect to one another, so that a particular background pixel/voxel, when not occluded by an organ, can be assumed to be of fixed intensity. Regarding the particular gated SPECT acquisition, the myocardium is the only object of concern. If the myocardium can be considered as a homogeneous object within a frame (i.e. constant in space), it is known that it becomes brighter from the end diastolic frame to the end systolic frame (i.e. variable in time) due to myocardial wall thickening and global uptake conservation [START_REF] Galt | Effects of myocardial wall thickness on SPECT quantification[END_REF]- [START_REF] Cooke | Determining the accuracy of calculating systolic wall thickening using a fast Fourier transform approximation: A simulation study based on canine and patient data[END_REF]. Therefore, we define one unknown constant myocardium value o n per frame in the following new balanced (it involves both the background and the object) criterion:

J 2 (c n , o n , b) = N n=1 Ω c 2 n (x) × {b(x) -S n (x)} 2 + Ω {c n (x) -1} 2 × {o n -S n (x)} 2 +β Ω ϕ b (|∇b(x)|) + γ N n=1 Ω ϕ c (|∇c n (x)|) • (8)

Constrained criterion: Active contour model

The algorithm developed in [START_REF] Kornprobst | Image sequence analysis via partial differential equations[END_REF] uses an initial estimate of c n (.) and makes it evolve iteratively pixel by pixel (the ϕ c -regularization term only provides a local regularization constraint, not a real geometrical constraint). Our aim is to propose a new method which (i) requires an initial segmentation for each frame (composed of one or several contours), (ii) interprets the inside of the contours as the regions of the moving (and warping) object, and (iii) makes the contours in each frame iteratively evolve toward the true moving object boundary by taking into account all of the frames. Therefore, we need to define a geometrically constrained criterion in which the segmentation unknowns are the moving object regions in each frame. We introduce one unknown object domain Ω n per frame. Let us denote: Ω n = Ω \ Ω n . The conditions of minimization of criterion (1) (see section 3) are still valid for criterion [START_REF] Kornprobst | Image sequence analysis via partial differential equations[END_REF]. Indeed: (i) If b(x) = S n (x) and c n (x) = 1, the two integrals of the data consistency term of (8) equal zero. This case corresponds to the background since b(x) = S n (x) (using the object domain notation:

x belongs to Ω n ). (ii) If b(x) = S n (x)
and c n (x) = 0, the first integral of the data consistency term equals zero, and the second integral equals o n = S n (x), which means that x belongs to the object (using the object domain notation: x belongs to Ω n ). We rewrite criterion (8) using the fact that Ω = Ω n ∪Ω n and Ω n ∩ Ω n = ∅:

J 2 (c n , o n , b) = N n=1 Ωn c 2 n (x) × {b(x) -S n (x)} 2 A + Ωn c 2 n (x) × {b(x) -S n (x)} 2 B + Ωn {c n (x) -1} 2 × {o n -S n (x)} 2 C + Ωn {c n (x) -1} 2 × {o n -S n (x)} 2 D +β Ω ϕ b (|∇b(x)|) +γ N n=1 Ω ϕ c (|∇c n (x)|) • (9)
We can delete the ϕ c -regularization term which ceases to be of use since c n (.) is no longer a function from Ω to the interval [0, 1] but a function from Ω to the set {0, 1}. Applying the above conclusions (c n (.) equals zero on Ω n and one on Ω n ), integrals A and D equal zero, thereby leading to the new criterion:

J 3 (Ω n , o n , b) = N n=1 Ωn {b(x) -S n (x)} 2 + Ωn {o n -S n (x)} 2 + β Ω ϕ b (|∇b(x)|) • (10)

Resolutions with respect to b(.) and o n

The minimization of criterion [START_REF] Yezzi | A statistical approach to snakes for bimodal and trimodal imagery[END_REF] with respect to b(.) and o n (Ω n being fixed) is achieved by solving the corresponding Euler-Lagrange equations. This criterion is a least square expression with respect to b(.) (if the ϕ bregularization term is omitted) and o n . Thus, it can easily be shown that the solutions are:

∀x ∈ Ω, b(x) = N n=1 ω n (x)S n (x) N n=1 ω n (x) if β = 0 (11) o n = S n (x), x ∈ Ω n ( 12
)
where ω n (x) equals one if x belongs to Ω n and is otherwise zero, and S denotes the average value of variable S. Equation ( 11) is a weighted time average. Equation ( 12) is a space average. If β is not equal to zero, the resolution with respect to b(.) is still a classical minimization problem [START_REF] Geman | Nonlinear image recovery with half-quadratic regularization[END_REF], [START_REF] Charbonnier | Deterministic edge-preserving regularization in computed imaging[END_REF].

Evolution velocity computation

Unregularized evolution velocity

The minimization of criterion [START_REF] Yezzi | A statistical approach to snakes for bimodal and trimodal imagery[END_REF] with respect to Ω n (b(.) and o n being fixed) is not straightforward. Indeed, the calculations of the Ω n Euler-Lagrange equations are delicate because the Ω n 's are integration domains. Moreover, these domains induce discontinuities: The intensity abruptly changes from o n to a value of background b(.) when crossing out of the boundary of Ω n . Then, computation of the derivative of ( 10) with respect to Ω n requires the use of the distribution theory to take into account the discontinuities. The idea is to embed Ω n into a family of evolving domains Ω n (t), t ≥ 0. Here, t represents the iteration index during computation of the solutions, which should not be confused with n, the frame index. Note that the ϕ bregularization term is independent of Ω n . Therefore, it is not involved in the following calculation.

Let us define:

F n (t) = Ω n (t) {b(x) -S n (x)} 2 + Ωn(t) {o n -S n (x)} 2 • (13)
Therefore:

J 3 (Ω n , o n , b) = N n=1 F n (t) + β Ω ϕ b (|∇b(x)|) • (14)
We would like to know how criterion [START_REF] Yezzi | A statistical approach to snakes for bimodal and trimodal imagery[END_REF] is modified if Ω n evolves slightly from instant t to instant t + dt. In other words, we would like to calculate

F ′ n (t) in the distributional sense. Let φ be a C ∞ -function with compact support. The distributional derivative of F n (t) is defined as < F ′ n (t), φ >. < F ′ n (t), φ(t) > △ = +∞ t=0 F ′ n (t)φ(t)dt = -< F n (t), φ ′ (t) > (property of the distributions) = - +∞ t=0 Ω n (t) {b(x) -S n (x)} 2 θ(t) φ ′ (t)dt - +∞ t=0 Ωn(t) {o n -S n (x)} 2 θ(t) φ ′ (t)dt • (15)
Integrating by parts:

+∞ t=0 θ(t)φ ′ (t)dt = - +∞ t=0 θ ′ (t)φ(t)dt+[θ(t)φ(t)] +∞ 0 • (16)
The last term of equation ( 16) equals zero since: φ(t = 0) = φ(t → +∞) = 0, due to the compact support of φ. Performing the same development with θ and by substitution into [START_REF] Siddiqi | Area and length minimizing flows for shape segmentation[END_REF], one gets:

< F ′ n (t), φ(t) > = +∞ t=0 θ′ (t)φ(t)dt + +∞ t=0 θ ′ (t)φ(t)dt = +∞ t=0 θ′ (t) + θ ′ (t) φ(t)dt • (17)
Equation ( 17) is true for any C ∞ -function φ with compact support. Then, we can conclude that: F ′ n (t) = θ′ (t) + θ ′ (t), in the distributional sense. Computation of θ′ (t) is not straightforward. However, it is a classical development in continuous media mechanics [START_REF] Haug | Methods of engineering mathematics[END_REF]. The final expression is:

θ′ (t) = Ω n (t) ∂ ∂t {b(x) -S n (x)} 2 + Γn(t) {b(x) -S n (x)} 2 v(x) • N (x)dΓ n (18)
where Γ n (t) is the oriented boundary of Ω n (t), v(.) is equal to (∂Γ n /∂t)(.) (therefore, it is the evolution velocity on Γ n (t)), and N (.) is the direct unit normal of Ω n (t).

Note: Active contours are naturally involved in the minimization of criterion [START_REF] Yezzi | A statistical approach to snakes for bimodal and trimodal imagery[END_REF] with respect to Ω n through the term v(x) involved during this development.

The integral term on Ω n (t) equals zero because {b(x)-S n (x)} 2 is independent of t. Similarly, since {o n -S n (x)} 2 is independent of t, θ ′ (t) equals:

θ ′ (t) = Γ ⋆ n (t) {o n -S n (x)} 2 v ⋆ (x) • N ⋆ (x)dΓ ⋆ n ( 19
)
where Γ ⋆ n (t) is the oriented boundary of Ω n (t), v ⋆ (.) = (∂Γ ⋆ n /∂t)(.) is the evolution velocity on Γ ⋆ n (t), and N ⋆ (.) is the direct unit normal of Ω n (t). Since Γ n (t) and Γ ⋆ n (t) represent the same contour with opposite orientations, one has: v(.) = v ⋆ (.), N (.) = -N ⋆ (.), and then:

θ ′ (t) = - Γn(t) {o n -S n (x)} 2 v(x) • N (x)dΓ n • (20)
Finally:

F ′ n (t) = Γn(t) {b(x) -S n (x)} 2 -{o n -S n (x)} 2 v(x) • N (x)dΓ n • (21)
The term v(x) • N (x) is the normal component of the evolution velocity, that is to say the useful component of the evolution velocity. To minimize criterion [START_REF] Yezzi | A statistical approach to snakes for bimodal and trimodal imagery[END_REF] with respect to Ω n , we have to find the expression of v(x)• N (x) which makes F ′ n (t) the most negative. However, we know neither the absolute value nor the sign of {b(x) -S n (x)} 2 -{o n -S n (x)} 2 . Let us denote this term A. The only way to guarantee that

F ′ n (t) is negative is to impose v(x) • N (x) to be equal to -A (in which case F ′ n (t) is the integral of -A 2 on Γ n (t)). Therefore: v(x) • N (x) = -{b(x) -S n (x)} 2 -{o n -S n (x)} 2 • (22)
As a conclusion, the evolution velocity on the boundary of Ω n is:

ν n (x) = ν n (x) N (x) (23) = {o n -S n (x)} 2 -{b(x) -S n (x)} 2 N (x) (24)
and the partial differential equation that models the active contour process is:

∂Γ n ∂t = ν n N • (25)
Assuming that b(.) and the o n 's are known, equations ( 24) and ( 25) can be interpreted as follows: If the contour is inside the object, only the second term of equation ( 24) contributes to the evolution since o n -S n (x) equals zero. The velocity is negative, therefore, outward contour motion occurs. Vice versa, if the contour is in the background (b(x) -S n (x) equals zero), only the first term contributes and inward contour motion occurs.

Note that from the beginning of section 4 to equation [START_REF] Galt | Effects of myocardial wall thickness on SPECT quantification[END_REF], the development is valid both for the 2-D case and the 3-D case:

• The sequence is either composed of 2-D frames or 3-D frames, and

• a contour is either a closed curve or a closed surface.

Curvature based regularization

The contour evolution process is driven by velocity ν n (.). The noise in the frames of the sequence can corrupt ν n (.). As a consequence, the segmentation contours might not be smooth and small parts of them may be left out during the evolution, thereby encircling noise peaks. A minimal length penalty is adapted to regularize such an evolution process [START_REF] Yezzi | A statistical approach to snakes for bimodal and trimodal imagery[END_REF], [START_REF] Mumford | Boundary detection by minimizing functionals[END_REF] since small sacrifices regarding the data consistency term minimization will allow smoothing of the main contours and shrinking of the small contours until they vanish (see sections 6.3 and 6.4 for more details). In the 2-D case, this penalty is expressed through the following criterion:

J 4,n (Ω n ) = Γn dp ( 26
)
where p is the arc-length parameter of curve Γ n . It can be shown [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF] that the minimization of ( 26) with respect to Γ n is achieved by solving the Euclidean curve shortening equation, or geometric heat equation:

∂Γ n ∂t = κ n N ( 27
)
where κ n (.) is the curvature of Γ n . This regularization term is added to criterion [START_REF] Yezzi | A statistical approach to snakes for bimodal and trimodal imagery[END_REF] after multiplication by a positive parameter δ to tune the relative weights of the two parts. Our variational criterion is then:

J 5 (Ω n , o n , b) = J 3 (Ω n , o n , b) + δ N n=1 J 4,n (Ω n ) • (28)
This regularization technique is also used in [START_REF] Yezzi | A statistical approach to snakes for bimodal and trimodal imagery[END_REF]. However, criterion (28) cannot be interpreted in the weighted length minimization framework (see section 6.4 for more details). Weighted length and area minimization flow derivations from first principles can be found in [START_REF] Yezzi | A geometric snake model for segmentation of medical imagery[END_REF], [START_REF] Caselles | Geodesic active contours[END_REF], [START_REF] Siddiqi | Area and length minimizing flows for shape segmentation[END_REF], and [START_REF] Angenent | On the Laplace-Beltrami operator and brain surface flattening[END_REF]. Therefore, the final partial differential equation to solve is:

∂Γ n ∂t = {ν n + δκ n } N • (29)
Extension of this curvature based regularization to the 3-D case is classical. Only the following changes are needed:

• Curve is replaced with surface,

• length is replaced with area,

• dp is replaced with dε, the infinitesimal area on Γ n , and

• curvature is replaced with mean curvature.

Level set implementation

The level set technique [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF]- [START_REF] Sethian | Level set methods[END_REF] greatly simplifies the management of the contour evolution (i.e. the resolution of equation ( 29)), especially for handling topological changes. The level set u n (.) is a function of x (2-D or 3-D point) with real values (u n : Ω -→ R) and is defined as follows:

• u n (.) is negative on Ω n ,

• u n (.) is positive on Ω n , and

• u n (.) equals zero on Γ n , i.e. for all p in [0, 1]:

u n (Γ n (p) x ) = 0 • (30)
The knowledge of the sign of u n (.) is equivalent to the knowledge of Ω n . Thus, the Ω n 's can be replaced with the u n (.)'s in criterion [START_REF] Haug | Methods of engineering mathematics[END_REF]:

J 5 (Ω n , o n , b) becomes J 5 (u n , o n , b
). In section 4.4.1, the Ω n 's are embedded into families of evolving domains Ω n (t). As a consequence, the u n (.)'s are also embedded into families of evolving level sets u n (., t). Equation ( 30) becomes:

∀x ∈ Ω, u n (x, t) = 0 • (31)
Differentiating [START_REF] Adalsteinsson | The fast construction of extension velocities in level set methods[END_REF] with respect to t (see the appendix), we obtain the following partial differential equation:

∂u n ∂t = {ν n + δκ n }|∇u n | (32)
where:

∀x ∈ Ω, κ n (x) = div ∇u n (x) |∇u n (x)| • (33)

Space-time segmentation algorithm

The algorithm was implemented in the 3-D case. Since criterion (28) must be minimized with respect to b(.), o n , and Ω n (n in [1..N ]) we use an alternate resolution scheme. Note that b(.) and the o n 's are also embedded into families of evolving terms, b(., t) and o n (t) respectively. Each iteration is divided into two steps. First, fixing domains Ω n , the minimization is achieved for b(.) and the o n 's. Second, fixing b(.) and the o n 's, the minimization is achieved for domains Ω n . The second step is performed by applying the active contour evolution velocities ν n (.) to the boundaries of domains Ω n , or equivalently to level sets u n (.). The segmentation algorithm is presented in table 3.

Initialization of the algorithm

The algorithm inputs are: b(., t = 0), o n (t = 0), and u n (., t = 0) (n in [1..N ]). b(x, t = 0) equals zero for all x. o n (t = 0) equals zero for all n. u 1 (., t = 0) is such that the zero level is a truncated hollow ellipsoid whose interior wall is not too far from the left ventricle interior wall in the first frame, and similarly for the exterior wall. u n (x, t = 0) = u 1 (x, t = 0) for all x and all n in [2.

.N ].

The parameters of the algorithm are: β, ϕ b , δ, and dt. We decided to set β to zero since the background of our SPECT sequences was not noisy enough to justify the use of the ϕ b -regularization. Then, the choice of ϕ b is not relevant. δ is chosen heuristically depending on the desired smoothness of the computed segmentation. If the segmentation is too irregular, δ must be increased. However, increasing δ too much may reduce the geometrical accuracy (see sections 6.4 and 6.5). Finally, dt is chosen such that the absolute value of the evolution velocity ν n (.) multiplied by dt is bounded by 0.1 and 0.2 during the first two or three iterations. These boundaries are arbitrary. They were chosen in order to guarantee the process stability.

Experimental evaluation

Our 3-D space-time segmentation method was tested on simulated and clinical 3-D gated SPECT sequences.

MCAT simulation

The method was tested on an MCAT simulation [START_REF] Pretorius | Evaluation of right and left ventricular volume and ejection fraction using a mathematical cardiac torso phantom for gated pool SPECT[END_REF], [START_REF] Lacroix | Evaluation of an attenuation compensation method with respect to defect detection in Tc-99m-Sestamibi myocardial SPECT images[END_REF]. Sixteen frames over the cardiac cycle were computed, but only height of them were used (every other frame starting from the first one). The frames were corrupted by a Poisson noise, thereby inducing signal to noise ratios in dB of: ∼32 in peak, ∼11.5 in energy, and ∼9.5 in variance. The myocardium region was not reformatted into a particular orientation. In the following, transverse slices refer to slices orthogonal to the body axis. The 128 3 generated volumes were cropped resulting in 45x40x40 images with a voxel size of 2.3x2.3x2.3 mm 3 . Only the left ventricle has been assigned a significant perfusion value, which is identical in each frame. Thus, there was no myocardium brightening effect during systole. The contrast of the myocardium with respect to the background was approximately 70/10. δ was 0.01 and dt was 4e-05. The algorithm was stopped after 300 iterations. The left ventricle segmentation is presented in figures 1, 2, and 3. It seems that systole extends from the first frame to the third or fourth frame with a rather regular contraction. Diastole seems to extend from the fourth or fifth frame to the eighth frame with a fast dilation in the first two diastolic frames, followed by a slower dilation in the subsequent frames.

For each frame, we have computed two segmentation error measures:

The number of voxels of the segmentation which does not belong to the MCAT simulated myocardium (oversegmentation), and the number of voxels of the MCAT simulated myocardium which are not included in the segmentation (under-segmentation). The mean oversegmentation error over the eight frames was 9.9% with a minimum and a maximum error of 7.4% (in the third frame) and 12.6% (in the eighth frame). The mean under-segmentation error over the eight frames was 6.3% with a minimum and a maximum error of 6.1% (in the fifth frame) and 6.8% (in the second frame).

Gated SPECT acquisition

The method was tested on standard clinical patient data obtained at rest using 99m Tc-Sestamibi with a triple-detector Picker Prism 3000 gamma camera (Picker International, Inc., Cleveland, Ohio, USA) equipped with low-energy, high-resolution fan-beam collimators. These data were acquired at the Centre Antoine Lacassagne, Nice, France. Sixty projections were acquired over 360 degrees using a frame rate of height frames per cardiac cycle. The height reconstructions were performed using a standard clinical filtered backprojection algorithm (Hamming filter). They were reformatted into short axis slices resulting in 64x64x23 images with a voxel size of 2.5x2.5x5 mm 3 . The maximum diameter of the left ventricle in the short axis slices was approximately 28 voxels. The contrast of the myocardium with respect to the background was approximately 45/12 in end diastole and 55/12 in end systole. δ was 0.015 and dt was 5e-05. The algorithm was stopped after 200 iterations. Only the left ventricle was segmented. The results are presented in figures 4, 5 and 6. It seems that systole extends from the first frame to the third or fourth frame with a rather regular contraction. Diastole seems to extend from the fourth or fifth frame to the eighth frame with a fast dilation in the first two diastolic frames, followed by a slower dilation which is almost imperceptible. In order to visualize the myocardium brightening effect, we have plotted the computed o n values against the frame number (see figure 7). This brightening effect is theoretically inversely proportional to the myocardium volume variation. This is the global uptake conservation principle. However, in practice, the value of the myocardium also depends on the resolution of the reconstruction, and then of the detector. The order of magnitude of the resolution of a SPECT detector is 1 cm. This is also the standard thickness of the heart in end diastole. This means that the value of the myocardium is certainly decreased in diastolic frames (increasing the brightening effect due to global uptake conservation). We have computed the products between o n and the number of voxels of the segmented myocardium in frame n. Using the level set, this number is the number of negative voxels of u n (.). The n products are theoretically proportional to the global uptake. They should be equal. The values obtained for the different frames had a standard deviation of 7.3% from the mean value with a minimum value in the eighth frame (end diastole) and a maximum value in the third frame (end systole). These results are consistent with the previous remark about detector resolution.

Ejection fraction computation

We have computed the ejection fractions of the MCAT simulation and the patient data. The method of calculation takes advantage of the level sets. At convergence, u n (.) is negative within the myocardium and positive everywhere else. Thus, the volume of the left ventricle is proportional to the number of positive voxels bounded by the negative voxels (the myocardium voxels) and a valve surface at the base of the myocardium. As suggested by [START_REF] Faber | Left ventricular function and perfusion from gated SPECT perfusion images: An integrated method[END_REF], we chose to model this valve surface with two planes. We have manually lo- where min and max are respectively the minimum and maximum number of left ventricle voxels over the sequence.

The ejection fraction of the MCAT simulation was 61.4% with a maximum left ventricle volume in the first frame and a minimum volume in the third frame. For the segmented myocardium, the maximum and minimum left ventricle volumes also occurred in the first frame and the third frame respectively. The computed ejection fraction was 60.3% corresponding to a 1.1% under-estimation (1.7% relative error).

We have also performed this computation on the healthy patient sequence of section 5.2. The maximum and minimum left ventricle volumes of the sequence segmentation occurred in the first frame and the fourth frame respectively. The computed ejection fraction was 57.5%, which is a normal value.

We have also segmented a gated SPECT sequence of a patient suffering from diffuse hypokinesis (results not presented in the paper). The ejection fraction obtained using ultrasound imaging was 30%. The maximum and minimum left ventricle volumes of the sequence segmentation occurred in the first frame and the fourth frame respectively. The computed ejection fraction was 27% corresponding to a 3% under-estimation (10% relative error).

The ejection fraction computation is an example of a dynamic feature measurement using a simple property of level sets: Their sign. See [START_REF] Zeng | Segmentation and measurement of the cortex from 3-D MR images using coupled-surfaces propagation[END_REF] for examples of static feature measurements (volume, shape index and thickness in-between two coupled level sets) on cortex segmentation using the sign and the value of level sets.

Discussion

Validity of the assumptions

We have made two assumptions. First, the background is assumed to be still. In 99m Tc-Sestamibi SPECT examinations, besides the myocardium, the liver, the kidneys (which are out of the field of view) and the spleen may have significant perfusion levels. These organs can be considered motionless if the motion due to breathing is not considered. The rest of the background is dominated by noise. As a consequence, this assumption is realistic in first approximation.

Second, within a frame, the myocardium value is assumed to be homogeneous. Before discussing this assumption, let us see if our results are in accordance. This is not necessary for the MCAT simulation. For the patient data, we have computed an intensity histogram of the segmented myocardium region for each frame of the sequence. These histograms are presented in figure 8. The homogeneity is not strict. This is mainly due to the smoothing of the reconstruction (near the edges of the myocardium) caused by the limited resolution of the detector. For each frame, we have also computed the number of voxels of the segmented myocardium whose values are in the interval [average value -15%,average value + 15%]. Table 4 presents the results. Between 45% and 57% of the voxel val- ues belong to this interval depending on the frame. As a conclusion, the assumption is fairly well respected. Nevertheless, this sequence is a healthy patient case for which the myocardium perfusion is correct. Figure 9 presents four situations of a myocardium with a hypoperfused region and the corresponding segmentation contours which would be obtained. In the first case (figure 9.a), the assumption is no longer valid. However, this case never appears in practice. In the second case (figure 9.b), the segmentation contours do not represent the myocardium itself. Nevertheless, the computed ejection fraction would still be correct. In the third case (figure 9.c), the computed ejection fraction would not be correct. It is necessary to perform tests using simulations and clinical sequences to determine the degree to which the ejection fraction deteriorates in such cases. In the fourth case (figure 9.d), the ejection fraction cannot be computed.

Advantages of the method

From a theoretical point of view, we have claimed that our method has the advantage of taking into account the entire sequence in each simultaneous frame segmentation process. Observing the active contour evolution velocity (equation ( 24)), this is not necessarily obvious. Indeed, only frame S n (.) is involved in the calculation of the evolution velocity ν n (.). Actually, all of the frames of the sequence are taken into account due to the presence of background b(.). At each iteration, b(.) is explicitly solved using all of the frames.

Generally speaking, what is the advantage of segmenting a sequence using a global process like the one we propose over several independent segmentation processes (one per frame)? To segment is to find homogeneous regions or boundaries (i.e. high gradients) between homogeneous regions. The previous question is equivalent to considering if it is easier (more efficient) to find high temporal gradients than high space gradients. The answer depends on the characteristics of the sequences and perhaps both types of gradient could be used in conjunction for maximal efficiency. Note that a temporal gradient and an optical flow are different approaches. An optical flow calculation assumes that a particular voxel x in frame n has the same value in frame n + 1 at location x + dx. Voxels are tracked from one frame to the next. A temporal gradient is the derivative of a voxel value with respect to time. The evolution velocity term {b(x) -S n (x)} 2 is the square of a temporal gradient first order approximation if b(.) is considered as a frame of sequence S n (.) at a time when the object was not yet in the field of view.

From a practical point of view, the use of the level set technique in our space-time segmentation process should greatly facilitate computation of dynamic features such as myocardial wall thickening. Indeed, at convergence, we have instantaneous information such as the myocardium region according to the sign of the level set. We have illustrated this by computing the ejection fractions of three gated SPECT sequences.

Experimental results

The segmentation errors (see section 5.1) on the MCAT simulation can probably be explained by one factor acting in two regions. This factor is the curvature penalty (see section 6.4). First, this penalty does not allow the active contours to fit highly concave shapes: This is the case in the left ventricle cavity, mainly in the apex region, in particular in the third and fourth frames. As a result, an over-segmentation occurs. Second, the MCAT myocardium includes voxels of the region of the valves. However, this region being too thin and having a lower perfusion level as well as being a highly convex region (which again cannot be accurately segmented because of the curvature penalty), it was not segmented. As a result, an under-segmentation occurs.

It seems difficult to overcome this problem. If there is noise, the contours may be irregular, which leads to segmentation errors. If we add the curvature penalty to smooth the contours (therefore, leading to fewer segmentation errors due to noise), the deep concavity or convexity problem arises possibly leading to segmentation errors.

These segmentation errors are presented in figures 10 and 11 for two values of delta (0.01 and 1000) and two levels of noise (levels 1 and 6, see table 5 in section 6.6). It can be observed that the over-segmentation errors are mostly located along the interior wall (in the apex region in particular) and that the under-segmentation errors are mostly located in the basal region.

Errors on the ejection fraction computations (section 5.3) may be due to two reasons. First, the quality of the segmentation. Second, the location of the valve planes. Indeed, we have located them manually. These computations were performed only to illustrate the possibility of obtaining, at reasonable cost, some dynamic features from the segmentation. Therefore, we did not make particular efforts to implement an ac-Figure 11: MCAT segmentation errors in the third frame (end systole). First and second row: Same configuration as in figure 10. curate and robust method for the valve plane location as proposed in [START_REF] Faber | Left ventricular function and perfusion from gated SPECT perfusion images: An integrated method[END_REF].

Algorithm convergence and curvature regularization

We have not studied the theoretical algorithm convergence. This is a delicate issue because of the alternate resolution scheme. Moreover, the unregularized evolution velocity ν n (.) can be positive on some parts of the contour and negative on some others. As far as we know, no work has been done yet to prove the convergence of an active contour process when the evolution velocity has not the same sign everywhere. Finally, note that the addition of the curvature based regularization does not allow us to interpret our method as a weighted area minimization (thereby preventing us from using directly the related theoretical results). Indeed, our criterion can be written:

n Ω A n + Γn dp , but not with the weighted area form: n Γn B n dp. Our motivation for using this particular regularization term was to compute a smooth segmentation and to avoid, while converging, formation of small contours around noise peaks. However, this curvature term may compete with the unregularized velocity to which it is added (see equation [START_REF] Pretorius | Evaluation of right and left ventricular volume and ejection fraction using a mathematical cardiac torso phantom for gated pool SPECT[END_REF]). For example, if a part of the contour is initially inside an object and the object has a convex shape, then the unregularized velocity is negative (see equation ( 24)), whereas the curvature will necessarily be positive when the contour gets closer to the object boundary. As a consequence, the regularized criterion does not allow accurate segmentation of deep concavities or deep convexities (the maximum allowed concavities and convexities depend on parameter δ: If δ is high, the maximum is low and vice versa). This possible competition between the two terms of the evolution velocity is probably another delicate point when studying the theoretical algorithm convergence.

The numerical test of convergence depends on the evolution of the zero level between two iterations. This evolution is measured using the evolution velocity values in the area of the zero level. If the evolution is less than an arbitrary threshold, the algorithm is stopped. However, since it is very unlikely that the segmentation contours go exactly through voxel centers, these evolution velocity values are never equal to zero and stabilize not necessarily close to zero. As a consequence, we simply stop the algorithm after a given number of iterations. For the sequences presented here, the computation time was 20 to 30 minutes on a 240 MHz HP workstation.

Influence of parameter δ

The weight of the curvature penalty, tuned by parameter δ, may influence the convergence of the algorithm. We have performed a series of segmentations of the MCAT simulation increasing δ from 0.01 to 1000. For each result, we have computed the mean segmentation errors described in section 5.1, the left ventricle volumes in the first frame and the third frame, and the ejection fraction. These results are presented in figures 12 and 13. As δ increases, the curvature penalty gets stronger and the segmentation contours converge toward closed surfaces of a smaller volume. Therefore, compared to the MCAT myocardium, the number of over-segmented voxels decreases and the number of under-segmented voxels increases. Moreover, the number of under-segmented voxels increases much faster (by 7.8% in all) than the number of oversegmented voxels decreases (by 2.0% in all). We have seen in section 6.3 that the region of concern for undersegmentation was the basal one, whereas the region of concern for over-segmentation was the interior wall, mostly the apex region. Since the curvature in the basal region is higher than the curvature in the apex region, the error depending on quality in the basal region (the under-segmentation error) is more dramatically influenced. However, until δ equals 50, both the oversegmentation error and the under-segmentation error remain quite constant (the over-segmentation belongs to [9.9%, 10.1%] and the under-segmentation belongs to [6.3%, 6.5%]).

Both first frame and third frame left ventricle volumes decrease with δ. The curvature penalty increasingly prevents the accurate segmentation of the apex region of the interior wall and the basal region where the curvature is even larger. This phenomenon in the basal region is more pronounced for the first frame than for the third frame since the myocardial wall is thinner in the first frame (so the curvature in the basal region is higher). As a consequence, the determined valve surface gets closer to the apex while δ is increased, the computed cavity volume then gets smaller, and the evolution is faster for the first frame. Indeed, in figure 13, it can be observed that the third frame left ventricle volume is quite constant until δ equals 300, whereas the first frame left ventricle volume decreases immediately after δ equals 100. As a consequence, the computed ejection fraction is reliable until δ equals 100 (absolute variation of 0.3% for δ between 0.01 and 100).

One can imagine that for a certain large δ, the data consistency term will be weaker than the curvature penalty. As a result, the segmentation contours will go through the actual myocardium boundary. This effect will occur first in the apex region where the myocardium is thinner. Moreover, it will occur first for the diastolic frames where the apical thinning is more prominent. This can be observed in figure 14 which represents the segmentation obtained when δ equals 2000. In addition to this phenomenon, note that the segmented myocardium is shorter (basal region) than when there is a much smaller δ, as is the case on figure 1.

Influence of noise

The noise level in the sequence may influence the convergence of the algorithm. We have performed a series of segmentations of the MCAT simulation, increasing the level of Poisson noise (see table 5 for the corresponding signal to noise ratios), while keeping δ equal to 20. In order to keep the convergence process stable while noise was increased, we introduced a smoothing step of the evolution velocity ν n (.). More precisely, after computation of ν n (.) and before its application to the level set u n (.), a three by three average filtering was performed on ν n (.). We prefer this heuristic technique to the one that consists of presmoothing the sequence (which is theoretically more justifiable) because we do not want to introduce a myocardium brightening effect due to the difference of myocardium thickness between systolic and diastolic frames. For each result, we have computed the mean segmentation errors described in section 5.1, the left ventricle volumes in the first frame and the third frame, and the ejection fraction. These results are presented in figures 15 and 16. The evolution of the segmentation errors when noise is increased is similar to what was observed for the influence of δ. In our simulations, it appears that the Poisson noise nibbles little by little the edges of the myocardium and the regions of low value such as the basal region. As a consequence, the noisy myocardium becomes smaller than the original one. Then, the segmentation contours converge toward closed surfaces of a smaller volume. The nibbling effect is more prominent in the basal region, which explains why the under-segmentation increases by 2.7% while the over-segmentation only drops by 0.9%. The first frame left ventricle volume decreases by 5.1% when noise is increased from level 1 to level 6. This can be explained by the nibbling effect in the basal region which causes the determined valve surface to be closer to the apex, reducing the computed cavity volume. The third frame left ventricle volume increases by 2.1% when noise is increased from level 1 to level 6. This is confusing. However, note that the nibbling effect at the basal region is less prominent in this frame since the myocardial wall is thicker than it is in the first frame. The nibbling effect on the edges of the cavity is certainly more prominent. As a consequence, the cavity volume slightly increases. This can be (barely) observed in figure 17 which presents the same slice of the segmentation obtained with noise levels 1 and 6. Going back to the first frame case, the nibbling effect on the edges should have the same consequence. However, the effect in the basal region is probably more pronounced. As a consequence, the ejection fraction keeps decreasing while noise is increased. The absolute variation is 0.8% between levels 1 and 2, and 3% between levels 1 and 6.

Figure 18 presents the segmentations of the first frame and the third frame obtained for levels 1 and 6. Note the small shortening of the segmented my- We have chosen to smooth the evolution velocity at each iteration. This is not computationally efficient, and should we have a proof of convergence of the algorithm, it would no longer be valid. However, this technique has the advantage (over the presmoothing of the sequence) of preserving the contrast of the sequence. Whereas presmoothing would probably decrease the contrast in diastolic frames where the myocardium has a lower value and is thinner. As a consequence, the segmentation task would be more difficult.

Finally, let us remark that the manner in which noisy simulations were obtained is not realistic enough. Monte Carlo simulations would be more effective. Nevertheless, the effects and their consequences seem acceptable.

Computation of b(.) and o n

Among the algorithm outputs, only the Ω n 's are relevant for the segmentation. However, b(.) and the o n 's must be computed since the evolution velocity applied to the Ω n 's depends on them. They appear to be auxiliary variables. The background is indeed not interesting in SPECT images. It is not taken into account for the diagnosis. We could make the same conclusion with respect to the o n 's. Using 99m Tc-Sestamibi, it is not possible to obtain absolute quantitative information. Nevertheless, the relative quantitative knowledge given by the o n 's may be used to compute an average myocardial wall thickening as suggested by [START_REF] Mochizuki | Assessment of systolic thickening with thallium-201 ECG-gated single-photon emission computed tomography: A parameter for local left ventricular function[END_REF] and [START_REF] Cooke | Determining the accuracy of calculating systolic wall thickening using a fast Fourier transform approximation: A simulation study based on canine and patient data[END_REF].

Future works

The assumption of a homogeneous myocardium value within a frame is reasonable. Hypoperfused regions will be excluded from the computed segmentation. Since SPECT is a functional modality, this is exactly what we are expecting. However, it is also desirable to segment the whole myocardium, whether it is homogeneously perfused or not. Thus, for example, hypoperfused region volumes could be computed. This anatomical segmentation using functional data requires the use of an anatomical model (like in [START_REF] Germano | Automatic quantification of ejection fraction from gated myocardial perfusion SPECT[END_REF], [START_REF] Nelson | A robust edge detection method for gated radionuclide ventriculograms[END_REF], and [START_REF] Faber | Left ventricular function and perfusion from gated SPECT perfusion images: An integrated method[END_REF]). We will try to incorporate such a model in the method using a constraint on the level sets.

Further experimental evaluation

The myocardium uptake varies from one patient to another due to different physiological properties. As a consequence, the contrast changes between clinical sequences. In a rigorous experimental evaluation context we should perform a sequence normalization prior to running the segmentation algorithm itself. This could be done by a linear mapping of the minimum to maximum voxel value range to the [0, 100] interval. Thus, the method would be independent of the original sequence contrast (particularly concerning the choice of dt).

The influence of the initialization of b(., t = 0), o n (t = 0), and u n (., t = 0) (n in [1..N ]) should be studied.

The influence of parameter δ has been studied in the MCAT simulation case. A more complete study including patient data should be done. b(.) and o n are obtained by computations of average values (a time average value for b(.) and a space average value for o n , see section 4.3). Therefore, the influence of noise on their solutions should be minimal. However, the influence of noise critically appears in the evolution velocity ν n (.). A preliminary study of the influence of noise has been done. It needs to be continued using Monte Carlo simulations and patient data.

The derivative of [START_REF] Adalsteinsson | The fast construction of extension velocities in level set methods[END_REF] κ n (.) is the curvature/mean curvature of Γ n , the zero level of u n (.). However, it is correct to extend its definition to any isolevel of u n (.). Then, for all x in Ω:

κ n (x) = div ∇u n (x) |∇u n (x)| • (40)
Similarly, ν n (.) is defined on Γ n . However, we need to compute it in a narrow band here and there from Γ n . We could use the evolution velocity extension proposed in [START_REF] Adalsteinsson | The fast construction of extension velocities in level set methods[END_REF] or [START_REF] Gomes | Reconciling distance functions and level sets[END_REF]. We chose to classically apply the formula obtained on Γ n within the narrow band [START_REF] Sethian | Level set methods[END_REF].

To solve equation [START_REF] Gomes | Reconciling distance functions and level sets[END_REF] we use a finite difference approximation of ∂un ∂t and a specific discretization scheme for |∇u n | (which is an unpublished evolution of the discretization scheme developed in [START_REF] Debreuve | Discretization schemes for volume segmentation using a level set method[END_REF]).

  (.) description of ĉn (x) below), b(.) is the unknown background image, ϕ b and ϕ c are regularization functions, and Ω is the frame domain. The minimization of (1) leads to an estimated background image b(.) and an estimated sequence segmentation ĉn (.). Without the ϕ bbackground regularization term (β = 0), b(.) reduces to a weighted time average of the frames S n (.): b

Table 3 :

 3 Segmentation algorithm b(., t = 0) ← initial background image o n (t = 0) ← initial myocardium value in frame n, n in [1..N ] u n (., t = 0) ← initial level set, n in [1..N ] Set parameters β, ϕ b , δ, and dt (see section 4.7) t ← 0 repeat b(., t) ← arg min b {J 5 with u n fixed to u n (., t), n in [1..N ]} (using (11)) for each n in [1..N ] do o n (t) ← arg min on {J 5 with u n fixed to u n (., t)} (using (12)) done for each n in [1..N ] do ν n (.) ← {o n (t) -S n (.)} 2 -{b(., t) -S n (.)} 2 κ n (.) ← div(∇u n (., t) / |∇u n (., t)|) u n (., t + dt) ← u n (., t) + dt {ν n (.) + δκ n (.)} |∇u n (., t)| done t ← t + dt t max ← t until convergence (see section 6.4) space-time segmentation ← zero level set of u n (., t max ),n in [1..N ]

Figure 1 :

 1 Figure 1: MCAT: 3-D space-time segmentation of the left ventricle. Chronological order from left to right and top to bottom.

Figure 2 :

 2 Figure 2: MCAT: Same transverse slice of the 3-D space-time segmentation. Chronological order from left to right and top to bottom.

Figure 3 :

 3 Figure 3: MCAT: Some transverse slices of the 3-D space-time segmentation of the 8th frame.

Figure 4 :

 4 Figure 4: Patient: 3-D space-time segmentation of the left ventricle. Chronological order from left to right and top to bottom.

Figure 5 :

 5 Figure 5: Patient: Base/apex median short axis slice of the 3-D space-time segmentation. Chronological order from left to right and top to bottom.

Figure 6 :

 6 Figure 6: Patient: Some short axis slices of the 3-D space-time segmentation of the 8th frame. From left to right and top to bottom: Base to apex.

Figure 7 :

 7 Figure 7: Patient: Computed o n values (in percent of the highest value) against the frame number.

Figure 8 :

 8 Figure 8: Histograms of the segmented myocardium region for each frame of the sequence. Chronological order is from left to right and top to bottom.

Figure 9 :

 9 Figure 9: Likely influence of a hypoperfused region (dark grey) on the segmentation (dotted line).

Figure 10 :

 10 Figure 10: MCAT segmentation errors in the first frame (end diastole). First row: A transverse slice with the over-segmented voxels in grey and the undersegmented voxels in black; From left to right, delta = 0.01 and noise level 1, delta = 1000 and noise level 1, delta = 20 and noise level 6. Second row: 3-D views showing both over-segmentation and undersegmentation errors; From left to right, same order as in the first row.

Figure 14 :

 14 Figure 14: MCAT: 3-D space-time segmentation obtained when δ equals 2000, first four frames.

Figure 12 :

 12 Figure 12: Influence of δ (horizontal axis) over the segmentation errors: Left, over-segmentation errors; Right, under-segmentation errors.

Figure 13 :

 13 Figure 13: Influence of δ (horizontal axis) over the ejection fraction: Left, left ventricle volume at end diastole (first frame); Middle, left ventricle volume at end systole (third frame); Right, ejection fraction.

Figure 15 :

 15 Figure 15: Influence of noise (horizontal axis) over the segmentation errors: Left, over-segmentation errors; Right, under-segmentation errors.

Figure 16 :

 16 Figure 16: Influence of noise (horizontal axis) over the ejection fraction: Left, left ventricle volume at end diastole (first frame); Middle, left ventricle volume at end systole (third frame); Right, ejection fraction.

Figure 17 :

 17 Figure 17: Influence of noise over the third frame left ventricle volume: A transverse slice of the segmentation obtained for noise levels 1 (left) and 6 (right).

Figure 18 :

 18 Figure 18: Influence of noise over the left ventricle segmentation. From left to right: First frame/level 1, first frame/level 6, third frame/level 1 and third frame/level 6.

Table 1 :

 1 Symbol list: Mathematical symbols and pa-

	rameters Symbol Category	Meaning
	R		The set of reals
	△ =		By definition, is equal to
	a • b		Dot product
	a, b		Duality pairing
	N	Input	Number of frames
			in a sequence
	α, β, γ	Parameters Constants of R +
	ϕ b , ϕ c	Parameter	Regularization functions
			from R + to R +
	n	Variable	Frame index in [1..N ]
	x	Variable	Pixel or voxel coordinates
	t	Variable	Algorithm iteration index
	dt	Parameter	Algorithm iteration step
	p	Variable	Arc-length parameter
	dp		Infinitesimal curve length
	dǫ		Infinitesimal surface area

Table 2 :

 2 Symbol list: Domains and unknowns

	Symbol	Category Meaning
	Ω	Input	Frame domain
	Ω × [1..N ] Input	N -frame sequence domain
	S n (.)	Input	Frame n of a sequence
	o n	Unknown Myocardium value in
			frame n (belongs to R + )
	b(.)	Unknown Background image of
			a sequence
	Ω n	Unknown Myocardium domain in
			frame n
	Γ n	Unknown Segmentation contour in
			frame n: Boundary of Ω n
	u n (.)	Unknown Level set of frame n whose
			zero level is Γ n
	N (.)	Variable	Direct unit normal of
			any isolevel of u n (.)
	ν n (.)	Variable	Evolution velocity
			magnitude on Γ n
	ν n (.)	Variable	ν n (.) N (.)
	κ n (.)	Variable	Curvature of any isolevel
			of u n

Table 4 :

 4 Number of voxels of the segmented my-

	ocardium whose values belong to [average value -
	15%,average value + 15%], in percent of the total num-
	ber of voxels of the segmented myocardium
	Frame	1	2	3	4
	Number of voxels 56.9 52.8 44.6 54.5
	Frame	5	6	7	8
	Number of voxels 53.8 55.7 56.1 55.8

Table 5 :

 5 MCAT: Signal to noise ratios in dB of the different levels of noise of the sequences (level 1 corresponds to the sequence presented in section 5.1)

	Noise level	1	2	3	4	5	6
	Peak	32.0 29.6 27.9 26.6 25.6 24.8
	Energy	11.5	8.3	6.6	5.4	4.4	3.6
	Variance	9.5	6.9	5.1	3.9	2.9	2.1

  with respect to t is: Using expression (29) of ∂Γ n /∂t, equation (36) is:∂u n ∂t = -{ν n + δκ n } N • ∇u n •(37)Since u n (.) is a level set and since it increases from the inside of Γ n to the outside of it, ∇u n (.) is normal to Γ n and is oriented toward the outside of Γ n . N (.) has been defined to have the opposite direction: N (.) = -∇un |∇un| (.). Then (37) is equivalent to:

						⇔	∂u n ∂t	= -	∂Γ n ∂t	• ∇u n •	(36)
						∂u n ∂t	= {ν n + δκ n }	∇u n |∇u n |	• ∇u n	(38)
					⇔	∂u n ∂t	= {ν n + δκ n }|∇u n | •	(39)
	∂Γ n ∂t	• ∇u n +	∂u n ∂t	= 0	(35)