Jérôme Declerck

Jacques Feldmar

Nicholas Ayache

Definition of a 4D continuous planispheric transformation for the tracking and the analysis of LV motion

Keywords: motion tracking, heart, LV, motion analysis, 4D

Cardiologists assume that analysis of the motion of the heart (especially the left ventricle) can give some information about the health of the myocardium. A 4D polar transformation is defined to describe the left ventricle (LV) motion and a method is presented to estimate it from sequences of 3D images. The transformation is defined in 3Dplanispheric coordinates (3PC) by a small number of parameters involved in a set of simple linear equations. It is continuous and regular in time and space, periodicity in time can be imposed. The local motion can be easily decomposed into a few canonical motions (radial motion, rotation around the long-axis, elevation). To recover the motion from original data, the 4D polar transformation is calculated using an adaptation of the Iterative Closest Point algorithm. We present the mathematical framework and a demonstration of its feasability on a series of gated SPECT sequences.

INTRODUCTION

Cardiologists assume that the analysis of the motion of the heart (especially the left ventricle) can give some information about the health of the myocardium. A huge effort has been made in medical image processing to track and analyse the motion of the LV, but due to the complexity of the modeling, this topic remains an open research problem.

Modern techniques provide 3D images which describe either the anatomy of the heart (MRI, US, for instance) or its functionality (Nuclear Medicine SPECT or PET imaging, for instance). It is possible to get sequences of such images over the whole cardiac cycle; such sequences are real 3D movies of the motion of the heart. The cardiac motion, like the motion of any real object must be therefore described as a 4D continuous and regular transformation of time and space. With the modality we are using for our experiments (gated SPECT), Corresponding author (e-mail: jdecler@sophia.inria.fr) the acquisition of an image lasts over several cardiac cycles; the gating procedure allows to reconstruct a 3D volume using information taken at similar times over different cardiac cycles, assuming that all the cycles are identical during the acquisition. In our protocol, we therefore suppose that the heart beats at a regular pace, this may not be the case for specific pathologies inducing irregularities of the cardiac pulse rate: these diseases are not supposed to be studied in this paper.

Many techniques have been proposed to track the LV motion. All of them attempt to find the correspondence between pairs of successive images. Most of the proposed methods in 3D define a model of the shape of LV surfaces (endocardium and/or epicardium), using classical snake-like models [START_REF] Mcinerney | A dynamic finite element surface model for segmentation and tracking in multidimensional medical images with application to cardiac 4D image analysis[END_REF][START_REF] Shi | Shape-based 4D left ventricular myocardial function analysis[END_REF][START_REF] Amini | Bending and stretching models for LV wall motion analysis from curves and surfaces[END_REF], spring-mass meshes [START_REF] Nastar | Vibration modes for non-rigid motion analysis in 3D images[END_REF] or more constrained generic surfaces such as freedeformed superquadrics [START_REF] Bardinet | Tracking and motion analysis of the left ventricle with deformable superquadrics[END_REF][START_REF] Bardinet | Tracking medical 3D data with a parametric deformable model[END_REF] or volumetric superquadrics [START_REF] Park | Analysis of left ventricular motion based on volumetric deformable models and MRI-SPAMM[END_REF][START_REF] Park | Deformable models with parameter functions: application to heart-wall modeling[END_REF]. The tracking is processed using conservation constraints based on proximity constraints [START_REF] Bardinet | Tracking and motion analysis of the left ventricle with deformable superquadrics[END_REF][START_REF] Bardinet | Tracking medical 3D data with a parametric deformable model[END_REF], differential properties of the surface [START_REF] Clarysse | 3D Boundary Extraction of the Left Ventricle by a Deformable Model with a priori Information[END_REF][START_REF] Shi | Shape-based 4D left ventricular myocardial function analysis[END_REF][START_REF] Amini | Bending and stretching models for LV wall motion analysis from curves and surfaces[END_REF][START_REF] Goldgof | Motion analysis of nonrigid surfaces[END_REF] or is directly computed from displacement or velocity information obtained in some specific MR imaging techniques: tags [START_REF] Radeva | Deformable B-Solids and Implicit Snakes for Localization and Tracking of MRI-SPAMM Data[END_REF][START_REF] Kraitchman | Semiautomatic tracking of myocardial motion in MR tagged images[END_REF][START_REF] Young | Tracking and Finite Element Analysis of Stripe Deformation in Magnetic Resonance Tagging[END_REF][START_REF] Denney | 3D displacement field reconstruction from planar tagged cardiac MR images[END_REF][START_REF] Park | Deformable models with parameter functions: application to heart-wall modeling[END_REF] or phase contrast [START_REF] Meyer | Tracking myocardial deformation using spatially-constrained velocities[END_REF][START_REF] Shi | A unified framework to assess myocardial function from 4D images[END_REF]. In other work, no shape model is computed: the tracking is processed directly from the volumetric image using optical flow methods [START_REF] Gorce | Estimation of three-dimensional cardiac velocity fields: assessment of a differential method and application to 3-D CT data[END_REF][START_REF] Song | Determining Cardiac Velocity Fields and Intraventricular Pressure Distribution from a Sequence of Ultrafast CT Cardiac Images[END_REF], conservation of differential elements of isophotes [START_REF] Benayoun | Dense non-rigid motion estimation in sequences of 3D images using differential constraints[END_REF] or using similarities of the intensity levels [START_REF] Thirion | Fast non-rigid matching of 3D medical images[END_REF].

Unfortunately, because the correspondence is defined between two successive images, regularity and periodicity in time is not guaranteed. Only a few studies impose temporal continuity or periodicity in their model: these studies deal with segmentation of 2D [START_REF] O'donnell | A periodic generalized cylinder model with local deformations for tracking closed contours exhibiting repeating motion[END_REF] or 3D images [START_REF] De Murcia | Reconstruction d'images cardiaques en tomographie d'émission monophotonique à l'aide de modèles spatiotemporels[END_REF][START_REF] Matheny | The use of three-and fourdimensional surface harmonics for rigid and nonrigid shape recovery and representation[END_REF][START_REF] Risler | A computer model for extracting moving heart surfaces from four-dimensional cardiac ultrasound images[END_REF]. Some other rare methods in 2D (Todd [START_REF] Constable | Development and evaluation of tracking algorithms for cardiac wall motion analysis using phase velocity MR imaging[END_REF][START_REF] Mceachen | A recursive filter for temporal analysis of cardiac motion[END_REF] or in 3D [START_REF] Thirion | Fast non-rigid matching of 3D medical images[END_REF][START_REF] Nastar | Vibration modes for non-rigid motion analysis in 3D images[END_REF] perform a posteriori time filtering. Moreover, all these tracking techniques ( [START_REF] Park | Analysis of left ventricular motion based on volumetric deformable models and MRI-SPAMM[END_REF]Park et al., 1994) excepted) do not provide intuitive parameters describing characteristic motions without non-trivial computation [START_REF] Bardinet | Tracking and motion analysis of the left ventricle with deformable superquadrics[END_REF][START_REF] Bardinet | Tracking medical 3D data with a parametric deformable model[END_REF][START_REF] Young | Three-dimensional left ventricular deformation in hypertrophic cardiomyopathy[END_REF]. The 4D polar transformation [START_REF] Declerck | Definition of a 4D continuous polar transformation for the tracking and the analysis of LV motion[END_REF] defined in this article aims to achieve four goals:

1. to define a class of transformation of time and space in which the temporal continuity and periodicity can be included, 2. to define a class of highly constrained transformations in order to have a relevant description of the LV motion with a reduced number of parameters, 3. to be able to retrieve canonical motions with minimal computation, providing an easy-to-interpret quantitative analysis of the motion.

4. last, but not least, to be a transformation which combines the unknown parameters in a linear way to make their estimation easier and robust.

We shall see that all these points are achieved in the mathematical formulation that we propose below. The paper is organised as follows: in section 1, we define the 4D polar transformation and the way to estimate it from a 4D (3D + time) data set. In section 2, a method is proposed to track in 4D the motion of the LV. Experiments have been conducted with a synthetic heart model and a gated SPECT sequence are presented in section 3. Section 4 draws conclusions concerning this work, its potential uses and future perspectives.

DEFINITION OF THE 4D PLANISPHERIC TRANSFORMATION

The idea of this study is to define a continuous and regular transformation of time and space. This transformation should also be adapted to describe with a minimum of parameters a complex motion such as the LV motion. This model of the deformation of the LV is a crude approximation compared to complex biomechanical models [START_REF] Hunter | The analysis of cardiac function: a continuum approach[END_REF] or highly descriptive kinematic models [START_REF] Arts | Description of the deformation of the left ventricle by a kinematic model[END_REF]. Given a point P (x, y, z) in cartesian coordinates and a time value t, the transformation gives a point Q (x 0 , y 0 , z 0 ) which is assumed to be the location of point P at time t. The cardiac motion is supposed to be regular in space and periodic in time. We therefore look for a differentiable function in spatial variables x, y and z and for a differentiable and a potentially periodic function in time variable t.

f : IR 3 IR ,! IR 3 P; t 7 ,! Q = f P; t
This definition of the 4D transformation yields the definition of 2 categories of functions which are easier to understand and which are intrinsically regular (figure 1): Definition of a 4D continuous planispheric transformation for the tracking and the analysis of LV motion 3

T P = f P;
In this paper, f is defined in order to describe locally some specific motions of points on the myocardium. We approximate the shape of the left ventricle as a stretched sphere in the long-axis direction. This is, of course, a very crude approximation as the shape of the heart is much more complex, however our goal is not a precise definition of the shape of the muscle, but a plausible discrimination of characteristic motions.

For that particular purpose, we separate the motion of a point of the heart into three canonical orthogonal motions (figure 2): motion 1: a radial motion which decribes the contraction or dilatation of the whole structure towards a "center" ; motion 2: a apico-basal rotation around the apicobasal axis which describes the twisting motion of the LV points ; motion 3: a motion (tangential to the surface r = C t ) which describes the elevation of the LV points in the apico-basal direction (the shortening along the long-axis during the systole). We describe these motions in a "3D-planispheric" coordinates (3PC) system, which is a combination of spherical and cylindrical coordinates. Our tranformation function is thus defined as a composition of three functions:

f = P2C F C2P
The function C2P switches from cartesian to 3PC, P2C switches back from 3PC to cartesian coordinates (of course, C2P = P2C ,1 ). F is the function which is described with the three basic motions in 3PC. The next two paragraphs detail the definitions of these functions.

Cylindrical or planispheric coordinates ?

The approach is inspired from [START_REF] Park | Analysis of left ventricular motion based on volumetric deformable models and MRI-SPAMM[END_REF][START_REF] Park | Deformable models with parameter functions: application to heart-wall modeling[END_REF]: in that study, the equations for the deformation of the superquadric model is expressed in cylindrical coordinates. In such coordinates, the decomposition of the local motion in our three canonical motions does not have the same relevance if the point where they are estimated belongs to a lateral wall (where the muscle is roughly cylindrical) or if the point is close to the apex (where the muscle is roughly spherical). For instance, a point belonging to a lateral wall and animated by an axial contraction does effectively contract towards the cylindrical axis (motion 1), but a point close to the apex animated with a similar motion does not contract, but undergoes a shift tangential to the wall. Thus, this motion is a twist (motion 2) or an elevation (motion 3) rather than a contraction (motion 1) (figure 3). We find it easier and more relevant to . P 1 belongs to a lateral wall, P 2 is close to the apex, both points are animated with motion 1. Solid arrows describe the motion in cylindrical coordinates (contraction towards the apico-basal axis) and dotted arrows describe the motion in spherical coordinates (contraction towards the center C). Motion 1 effectively describes a contraction for P 1 , but not for P 2 . use 3D-planispheric coordinates rather than cylindrical ones to decompose the local motion of points of the myocardium.

3D-Planispheric coordinates C2P

: in 3D cartesian space, we define a 3D-planispheric reference system given a center C, a base B and a set of two orthogonal vectors u and v (figure 7). In order to fit with our description of the heart, we choose u as a vector parallel to the apico-basal direction, and v parallel to the septo-lateral direction. C is chosen in the center of the cavity, and B in the center of the base.

For each point P (x, y, z), a center point H P is defined on line (CB). From this center point, a distance and two angles (latitude θ and longitude φ) are calculated just as in the classical spherical coordinate system. In the spherical system, H P is the center C. In the cylindrical system, H P is the orthogonal projection of P on the line CB. Our purpose here is to define a combination of both spherical and cylindrical coordinate systems (figure 4), in order to describe the position of P in "roughly" spherical coordinates around the apex (where the shape of the LV is roughly spherical) and in "roughly" cylindrical coordinates around the base (where the shape of the LV is roughly cylindrical). In our system, the position of H P on the line CB is given by the simple formula: For low values of θ (P around the apex), H P is close to C and shifts away from C with a distance increasing with θ2

CH P = 1 ,cosθCB (1)

2

: around the apex, the 3PC system is thus close to the spherical one. For θ around π 2 (P around the apex), H P is close to B, the distance BH P varies linearly with θ , π

Of course, ( 1) is an implicit formula: H P gives the angle θ, but we need θ to locate H P . appendix A details the method we have developed to compute the location of H P given a point P in space.

In our 3PC system, a surface of constant r is represented as a disk in a plane, like in a classical map projection in topography (this is why we use the word "planispheric"). The coordinates X, Y and R in this system are defined as follows:

X = θ π cosφ Y = θ π sinφ (2) R = r σ r
where σ r is a normalization coefficient so that X, Y and R are dimensionless and vary within a similar range of values. cartesian coordinates and the (X, Y, R) coordinates in the 3PC system. In this system, apex point (X = 0, Y = 0) is the "south pole" of the projection, the points on the circle X 2 + Y 2 = 1 (θ = π) are the same cartesian point, the "north pole" of the projection, featuring a point in the direction of the base or of the aorta. Around this point, the distorsion between the cartesian and our planispheric representation is maximum, but there should not be any cardiac points in this area (figure 6).

South pole

Lat.

Inf. myocardium North pole

Sep.

B C

Lat.

North pole circle

Ant.

Sep.

Figure 6.

Left, the myocardium. From the center C, the limit of the basis draws a cone (dark gray) around the north pole. Right, in the 3D-planispheric map, the left ventricle (light gray) appears like a plate, the cone is a circular stripe around the heart.

P2C : conversely, given a point (X, Y, R) in the 3PC system so that X 2 +Y 2 1, we can compute its cartesian coordinates (x, y, z) by calculating θ = π p X 2 + Y 2 and φ without ambiguity with the expressions of cosφ and sinφ. The center H P is calculated with (1).

The coordinate system we use is similar to the prolate spheroidal coordinates (PSC) system described in the litterature [START_REF] Waks | Cardiac motion simulator for tagged MRI[END_REF]. The equations are presented in appendix B, figure 7 shows ten surfaces of constant R in each of the two systems. There are three minor differences between the two systems:

1. the surfaces of constant R (for R values around what they should be to describe an average LV) are narrower around the apex than in the 3PC system: the shape of these surfaces is closer to the shape of an average LV;

2. there is an interval of R for which surfaces of constant R are close to a shape of an average LV. For those values, the location of point C is closer to the apex in the PSC system than in the 3PC, potentially yielding to forbidden intersection of segment [CB] with the myocardial wall;

3. the surfaces of constant θ are cones in the 3PC system, they are confocal one-sheet hyperboloids in the PSC system. The local coordinates are orthogonal in PSC, whereas they are approximately orthogonal in 3PC. On the other hand, the decomposition of the motion is less intuitive in the PSC system.

In the sequel of this article, we concentrate on the 3PC system in order to obtain a closer approximation of the LV shape and a more intuitive decomposition of the motion.

The function in 3D-planispheric coordinates

In the 3D-planispheric system, given a point

P (X, Y, R), the transformed point Q (X 0 , Y 0 , R 0 ) through F is expressed as follows: X 0 = a 0 X ,a 1 Y + a 2 Y 0 = a 1 X + a 0 Y + a 3 (3) R 0 = a 4 R + a 5
X 0 and Y 0 are defined by a similarity function applied to X and Y. A similarity is a combinate of a 2D rotation by angle α, a uniform scaling of ratio k and a translation. R 0 is defined as an affine function of R. The similarity and the affine parameters a p (p = 0 : : : 5) are continuous and differentiable functions of r, θ, φ and t.

Defining the transformation in the 3PC system allows to write: a) linear expressions in the parameters a p , b) a simple computation of the canonical motion decomposition (radial motion, rotation, elevation) from the a p and, c) a very compact description of the deformation.

Analysis of the motion in the canonical decomposition

Our canonical motions are retrieved with the following formulae:

1. the radial motion ratio (motion 1) is given by: R 0

R = a 4 + a 5 R (4)
2. the linear relationships between X 0 , Y 0 and X, Y define a 2D similarity such that

k = q a 2 0 + a 2 1 (5) α = atan2 a 1 k ; a 0 k (6)
α is the rotation around the long-axis (motion 2) in 3PC;

3. k is the scale factor corresponding to an elevation magnification in latitude, which is our motion 3. In our display, we compute θ 0 ,θ θ .

The 4D polar transformation is defined once the parameters a p are determined. Because of the simplicity of ( 5) and ( 6), it is possible to easily analyse the motion using the parameters a p . Because the variation of the parameters a p is smooth and regular with variables r, θ, φ and t, the parameters which describe our canonical motions are also smooth and regular in time and space.

Degrees of freedom of the parameters, time dependency as a hard constraint

In order to define a smooth and continuous 4D transformation, the parameters depend on the location of the point and the instant at which the transformation is calculated. In our formulation, we choose the parameters as polynomial functions in r and θ and quadratic periodic B-splines [START_REF] Risler | A computer model for extracting moving heart surfaces from four-dimensional cardiac ultrasound images[END_REF][START_REF] Farin | Curves and Surfaces for Computer Aided Geometric Design[END_REF] in φ and t:

a p r; θ; φ;t = n r ,1 ∑ i=0 n θ ∑ j = 0 j 6 = 1 n φ ,1 ∑ k=0 n t ,1 ∑ n=0 A p i; j;k;n 0 @ r σ r 1 A i 0 @ θ π 1 A j B Φ k φ B T n t (7)
for p = 0 : : : 4. If we keep for a 5 an expression like (7), a 4 and a 5 are correlated because R = r σ r . We therefore simplify a 5 as follows:

a 5 r; θ; φ;t = n θ ∑ j = 0 j 6 = 1 n φ ,1 ∑ k=0 n t ,1 ∑ n=0 A 5 j;k;n 0 @ θ π 1 A j B Φ k φ B T n t (8)
with the following notation:

n r is the number of parameters which define the polynomial function of variable r: the degree of this polynomial is n r ,1.

n θ is the number of parameters which define the polynomial function of variable θ: to be differentiable in points for which θ = 0, the polynomial must have no term in θ (a p θ = a 0 p + a 2 p θ 2 + a 3 p θ 3 : : : ). a p is therefore a polynomial of θ of degree n θ . n φ is the number of control points of the B-spline periodic curve of variable φ. B Φ are the B-spline basis functions associated to a classical regularly distributed 2π-periodic set of knots. n t is the number of control points of the B-spline curve of variable t. B T are the B-spline basis functions associated to a classical regularly distributed set of knots, this basis can be periodic or not.

The originality of the transformation is in the fact that the continuity and potentially the periodicity in time is a "hard" constraint. We can implicitly look for time-periodic transformations.

Using quadratic B-splines (with a set of regularly distributed knots in our current implementation) ensures C 1 continuity in φ and t; the function a p r; θ; φ; is a (potentially periodic) piecewise polynomial. Due to the definition of Bsplines, the influence of possible outliers remain local [START_REF] Risler | A computer model for extracting moving heart surfaces from four-dimensional cardiac ultrasound images[END_REF].

To ensure the continuity for θ = 0, we must impose the constraint A p i;0;k;n = A p i;0;0;n for each k. There are thus n φ , 1 equations for each i and n. Finally, we get a number of control points N CP = 5:n r + 1:n φ :n θ ,1+1:n t .

The transformation is completely defined given a center C and two orthogonal vectors (they define the 3PC system) and a set of N CP control points (real numbers) A p i; j;k;n .

Estimation of a 4D planispheric transformation 1.4.1. The least squares criterion

Having a set of matches (P l , Q n;l ) for different times t n (n = 0 : : : T ,1), we define a least squares criterion to estimate the 4D planispheric transform which could "best" fit the list of matches:

8n 2 0 : : : N ,1; 8m 2 0 : : : M ,1;

f P l ; t n ' Q n;l (9)
The least squares criterion is therefore written:

J f = T,1 ∑ n=0 N,1 ∑ l=0 α n;l : d f P l ; t n ; Q n;l 2 (10)
where d; is the distance and α n;l is the weight related to the reliability of the match (P l , Q n;l ).

If we choose the euclidean distance for d in cartesian coordinates, the criterion is not quadratic in the A p i; j;k;n , and its derivatives with respect to the A p i; j;k;n are very difficult to linearize. We prefer to choose for d the euclidean distance expressed in 3PC (X 2 +Y 2 + R 2 ). The criterion is then quadratic in the A p i; j;k;n . d is a distance if and only if X 2 +Y 2 1 (i.e. iff (X,Y) does not belong to the circle of the "north pole" of the planispheric map). As the center C is well inside the cavity and the base point B is in the center of the base circle, we are sure that all data points remain in a "security" cylinder in the 3PC system (maximum expected value for θ is around π=2, so X and Y are less than 0:5. See figure 6).

Minimization of the criterion

The criterion expressed with d as the euclidean distance in 3PC is quadratic in the control points A p i; j;k;n . Differentiating it with respect to A p i; j;k;n gives a linear system which is solved with a classical conjugate gradient method. The size of the matrix is 5:n r + 1:n φ :n θ , 1 + 1:n t . In fact, the linear system can be split in two independent subsystems, one for a 0 , a 1 , a 2 , a 3 and the other for a 4 and a 5 .

Assembling the matrix is an operation which obviously depends linearly with the N, the number of matches. Solving the linear system is an operation which depends on N CP , the size of the matrix : a classical conjugate gradient needs ON CP 2 elementary operations per iteration. In our experiments, the ON operation is much more costly than the ON CP 2 one (ratio is approximately 5:1).

TRACKING THE 4D MOTION OF THE LV

We define in this section an adaptation of the Iterative Closest Point algorithm [START_REF] Besl | A method for registration of 3D shapes[END_REF][START_REF] Zhang | Iterative point matching for registration of freeform curves and surfaces[END_REF] which gives an estimation of those matches for the least squares minimization: it is possible this way to calculate the "best" function with respect to a distance criterion.

The motion is tracked in a heart image sequence (in our experiments, gated SPECT). Points featuring the edges of the heart are extracted and matched. The result of the matches between pairs of points in the images of the sequence is used to estimate a 4D polar transformation.

Matching the feature points

The matching method is an enhancement of the iterative closest point [START_REF] Besl | A method for registration of 3D shapes[END_REF][START_REF] Zhang | Iterative point matching for registration of freeform curves and surfaces[END_REF][START_REF] Feldmar | Rigid, affine and locally affine registration of free-form surfaces[END_REF], adapted to our problem. Given a point P (x, y, z) in cartesian coordinates and a time value t, the transformation gives a point Q (x 0 , y 0 , z 0 ) which is assumed to be the location of point P at time t. To estimate a 4D planispheric transformation f , we therefore need to know the matches be-tween points P l of the first image (t = 0) and points Q n;l of the image at time n (t = t n ). We thus look for f so that

f P l ; t n ' Q n;l
We define a criterion:

C f = T,1 ∑ n=0 N,1 ∑ l=0 α n;l : d f P l ; t n ; CP n f P l ; t n 2 : ( 11
)
CP n calculates the closest point to a 3D point among the feature points FP n extracted in image n. We use kd-trees to compute this function [START_REF] Preparata | Computational geometry, an introduction[END_REF]. The criterion is the sum of the residual distances between the estimated location of the points P l at time t n and the feature points extracted in Image n. α n;l is a weighting coefficient which depends only on (P l , t n ) (figure 8). In our experiments, the model is the list of data points at end diastole. . The matching criterion measures the sum of distances between f (P l , t n ), the estimated point P l at time t n and Q n;l , the closest feature point of Image n, for all existing P l and t n .

Minimizing the criterion

The minimization process is iterative, given an initial transformation f 0 . This initial transformation is chosen in our experiments to be the identity (nothing moves anywhere !). Each iteration k splits into three steps:

1. For each 4D point (P l , t n ), we calculate f k,1 P l ; t n (which should be the location of point P l at time t n ) and we identify its closest feature point CP n;l in Image n. We therefore end up with a list of possible pairs of matched points.

2. For each time t n and for each type of boundary (endocardium, epicardium), we calculate the residual distance kf k,1 P l ; t n ,CP n;l k for each pair, and we decide whether a pair is reliable or not: we first eliminate pairs for which the residual distance exceeds a fixed threshold. Second, we compute the mean µ and the standard deviation σ attached to the remaining pairs. We then eliminate the points for which the distance is greater than another threshold depending on the distance distribution (µ + c:σ, where c can be easily set using a χ 2 table [START_REF] Feldmar | Rigid, affine and locally affine registration of free-form surfaces[END_REF][START_REF] Feldmar | Recalage rigide, non-rigide et projectif d'images médicales tridimensionnelles[END_REF]).

We get for this iteration a list S k of reliable pairs of matched points. Notice that if a point is not matched in this iteration, it may be matched in one that follows.

3. With the filtered list S k of pairs of points, we calculate f k which is the best least squares fit for the pairs of points.

The iterative process stops when a maximum number of iterations is reached, or when S k = S k,1 . [START_REF] Feldmar | Rigid, affine and locally affine registration of free-form surfaces[END_REF][START_REF] Feldmar | Recalage rigide, non-rigide et projectif d'images médicales tridimensionnelles[END_REF] gives further details about this adaptation of the iterative closest point algorithm, for instance about the convergence properties.

Definition of the closest point

The matching function CP n takes into account for each point its geometric position and the local direction of the intensity gradient calculated while extracting the edges. Considering 2 oriented points (P;n P ) and (Q;n Q ), where n P and n Q are the directions of the intensity gradient at point P and Q respectively, the distance between them is calculated as follows:

dP; Q 2 = α:kPQk 2 + kn P ,n Q k 2 (12)
where α is a weighting coefficient for normalisation. The local direction of the gradient defines which border a points belongs to: if the gradient is oriented towards the center, the point belongs to the endocardium. If not, it is assumed to belong to the epicardium. This separation avoids mismatches between points of two different boundaries and speeds up the computation: one kd-tree is more costly to manipulate than two kd-trees of half size. This double definition of a point (location + direction) refines the matching criterion and makes it more robust and precise. We know that such features must be used with caution, especially when trying to define a distance between two features [START_REF] Pennec | Randomness and Geometric Features in Computer Vision[END_REF]. However, the formula ( 12) must be written as a sum of squares in order to keep the convergence properties of the process [START_REF] Feldmar | Rigid, affine and locally affine registration of free-form surfaces[END_REF][START_REF] Feldmar | Recalage rigide, non-rigide et projectif d'images médicales tridimensionnelles[END_REF].

Computing an optimal 3D-planispheric coordinate system

A keypoint in the estimation of the 4D transformation is the definition of a 3PC system (a center, a base, an apico-basal vector, a septo-lateral vector and a normalisation factor σ r ).

The coordinate system

In [START_REF] Declerck | Automatic registration and alignment on a template of cardiac stress & rest SPECT images[END_REF], we define a method to align a SPECT heart image with a template using a non-rigid transformation. This method gives a transformation from the normalized coordinates of the template to the patient's case.

As the transformation which deforms the template is sufficiently free (B-spline tensor product), the template can be chosen as a rough approximation of a LV. Here, we choose two truncated ellipsoids (one for the endocardium and one for the epicardium). The parameters have been set manually not to design a precise shape: the idea is just to have a "goodlooking" one. We define for this template a center, a basis point and a point in the lateral wall so that all three define a reference system [C, u, v].

The template is matched with the edges of the image of the heart at end diastole (largest volume). With this transformation, we deform the reference system of the template to the patient case (figure 9). Calling S (for "shape") the spline transformation deforming the template to the patient's case, the reference system is defined as follows:

For the template C P L P , u P jC P L P : u P k::: k (so that u P :v P = 0) -infero-anterior vector : w P = u P v P

-center: C T -base: B T -lateral point : L T -apicobasal vector : u T = C T B T kC T B T k -septo-lateral vector : v T = C T L T kC T L T k (L T is such that u T :v T = 0) -infero-

Choosing σ r

The normalisation factor σ r is used to make the R coordinate dimensionless, as are the other coordinates X and Y. Changing σ r changes the shape of a surface in the planispheric geometry by a scaling in the R direction (the lower σ r , the "higher" the surface). The closest point in this surface to a given point P varies with σ r (figure 10):

When σ r approaches 0, the R value becomes very large compared to X and Y, Q 0 , the closest point to P tends to a point with the same R (figure 10, left). This implies that if we use those matches for the least square criterion on distances, the tangential motions (those which change only X and Y) are privileged and the radial motion (those which change only R) becomes negligible. When σ r tends to infinity, the R value becomes very small compared to X and Y. Q ∞ , The closest point to P is a point with same X and Y (figure 10, right). This implies that if we use those matches for the least square criterion on distances, the radial motion is privileged and the tangential motions are negligible.

Giving a value to this factor therefore amounts to choosing a weighting between purely tangential and purely radial motions.

As the latitude of the basal points approaches π=2, the absolute values of X and Y do not exceed 0:5. For an average heart, it appears that the maximum distance (meaning the r value) of a point of the myocardium to the axis does not exceed CB. We thus choose σ r = 2:CB, so that R does not exceed 0:5 as for X and Y.

o oo P o Q Q 0 R Y X S Y R Y S 0 R X Q X S Q

EXPERIMENTS

We present here experiments conducted on a series of gated SPECT image sequences provided by Pr. M.L. Goris, Stanford University Hospital (California, USA). There are 8 images in the sequence, the size of the images is 64x64x64, pixel size is 2.5mm isotropic. The temporal sampling is uniform and covers the entire cardiac cycle.

Extraction of feature points

Each image of the sequence is resampled in the polar geometry defined in [START_REF] Declerck | Automatic registration and alignment on a template of cardiac stress & rest SPECT images[END_REF]. This reference describes a method to extract edges in nuclear medicine myocardial perfusion images, we recall here the main ideas: in a polar geometry with a center well inside the cavity, the heart looks like a thick plate. We detect edges in this image with a Canny-Deriche recursive filter [START_REF] Monga | 3D edge detection using recursive filtering: application to scanner images[END_REF]. Cardiac boundary points are easily detected and filtered in this geometry: in an ideal situation, where the heart is a volume of pixels with high intensity values in an image with a high signalto-noise ratio, starting from the center of the cavity along a radius, the first edge is assumed to belong to the endocardium, the following edge is assumed to belong to the epicardium (figure 11). This heuristic constraint is used to filter out the non-cardiac edges (in SPECT images, there may be some parts of the liver, or artificial edges generated by binary masks). In [START_REF] Declerck | Automatic registration and alignment on a template of cardiac stress & rest SPECT images[END_REF], the method is tested on a database of 40 pairs of images, the proportion of erroneous edges was estimated satisfactory enough for the matching process to work. Figure 12 shows the result of the extraction of feature points on four different times of a gated SPECT sequence of a healthy heart.

Retrieval of trajectories

From these data points, we look for a periodic 4D polar transformation with n r = 2, n θ = 3, n φ = 6, n t = 5 (which gives 715 control points, for roughly 50000 data points). Figure 14 shows the surface of a healthy LV (VANN) deformed by the 4D polar transformation. This is the surface of a template matched with the method defined in [START_REF] Declerck | Automatic registration and alignment on a template of cardiac stress & rest SPECT images[END_REF]. Because the feature points we extract from the images are not structures, it is difficult to visualise them efficiently. The trajectories of the points are smooth and periodic, as illustrated by figure 13. Figure 15 shows the surface of a pathological LV (BOJO, suffering from septal akinesia due to an infarct) deformed by another 4D transformation: the akinesia is revealed by the relative short trajectories of points around the septum (on the left in each frame of figure 15).

Extracting feature points takes around 20 seconds per image on a DEC Alpha 400 MHz workstation. The computation of the 4D planispheric transformation takes around 15 minutes on the same machine.

Analysis of the motion

We show in figure 16 a display of the physical parameters values on the surfaces of the heart. On the top, the normal VANN case and on the bottom, the pathological BOJO case. For the radial motion, the values vary between 0 and 40 %, for the elevation, between -5 % and 5 % and for the apico-basal rotation, from -10 to 10 degrees. The values are computed at end systole.

In figure 17, we show a display of the amplitude and the phase of the first harmonic of the radial motion for three normal hearts (CRIR, HUVM and VANN) and a pathological one (BOJO). This crude Fourier analysis shows how the parameter varies over the whole sequence. For the normal cases, there is a variability in the computed values of the parameters:

1. around the base, the values are hasardous due to a bad identification of the base, potentially leading to mismatches. A better segmentation of the base should avoid this kind of behaviour;

2. apart from this region, the distribution of amplitude values is coherent for the three cases, showing that the septum contracts less than the other walls and that the contraction is more important on the endocardium than on the epicardium. The amplitude reveals potential akinetic areas;

3. the phase reveals potential dyskinetic regions. The values distribution shows the septum contracts slightly later than the other walls. This is particularly visible on cases VANN and HUVM.

However, those values must be interpreted with caution. As the matching procedure works with a closest point method, the tangential motions which are retrieved are not fully reliable (elevation and apico-basal rotation). As a matter of fact, only the radial motion (an "orthogonal" motion) matches what can be expected from a healthy heart [START_REF] Park | Analysis of left ventricular motion based on volumetric deformable models and MRI-SPAMM[END_REF][START_REF] Young | Three-dimensional left ventricular deformation in hypertrophic cardiomyopathy[END_REF] and figure 17 shows, for a pathological case, the sensitivity of this parameter: amplitude: the septum appears in dark blue, featuring low values. The akinesia is thus located and quantified; phase: around the septum, the values distribution shows strong variations, due to the fact that the amplitudes are low in this area: the approximation of a nearly flat curve (showing the evolution of the parameter over time) by a sinusoid is therefore very sensitive.

Perspectives for the validation

The motion we are able to retrieve and analyse using gated-SPECT images sequences must be interpreted with caution. The amplitude shows how wide the motion is, the phase shows the synchronisation.

Due to the low resolution of the images, it is difficult to get a precise information. Second, any tangential motion cannot be reliably retrieved using only feature-based techniques without markers. The parameters we are able to compute may be useful if there is a possibility to demonstrate that they can be used for a detection of a pathology, by separating normal and abnormal hearts into two statistically different classes. This validation should be processed on a dataset of heart images of which the pathology or healthy state is known. For a given database, the sensitivity and the specificity can be calculated and can show the usefulness of our approach on a quantitative basis.

Another way to validate our decomposition of the motion is to check that it corresponds to a real motion. Tagged MRI yields images in which the motion of soft tissues at a number of discrete points is easily detectable, can be measured [START_REF] Mcveigh | MRI of myocardial function: motion tracking techniques[END_REF][START_REF] Kraitchman | Semiautomatic tracking of myocardial motion in MR tagged images[END_REF][START_REF] Young | Tracking and Finite Element Analysis of Stripe Deformation in Magnetic Resonance Tagging[END_REF][START_REF] Young | Three-dimensional left ventricular deformation in hypertrophic cardiomyopathy[END_REF][START_REF] Denney | 3D displacement field reconstruction from planar tagged cardiac MR images[END_REF] and then compared to our computed motion.

These two validation processes are currently under study, partial results have been obtained [START_REF] Declerck | Étude de la dynamique cardiaque par analyse d'images tridimensionnelles[END_REF] and will be the subject of a forthcoming article.

CONCLUSION

In this work, the mathematical framework for a new class of transformation is defined: a 4D planispheric transformation is a differentiable function in space and time coordinates and potentially periodic in time. A small number of parameters constrain the definition of the function and there is a simple relationship between the estimated parameters a p and the "canonical" motions defined for a moving LV (radial motion, rotation, elevation). We demonstrated the feasability of the method on a series of gated SPECT sequences.

This will be the basis for a number of experimental studies both on nuclear medicine and tagged MR data in collaboration with Pr. Michael Goris (Stanford University Hospital) and Dr. Elliot McVeigh (Johns Hopkins University).

In order to refine the tracking procedure, we are also working on defining feature points inside the myocardium. Those points added to the edges we have already defined will give landmarks in the entire myocardial volume and not only on its boundary.

A. 3D PLANISPHERIC COORDINATES

This section is dedicated to the problem of finding H P given a point P in space. In order to avoid cumbersome notation, we rename H P as H. The problem is, given two points C and B and a point P, to find a point H on the line CB so that:

We look for λ 0 so that

f λ 0 = 0 (15)
To solve this equation, we use a Newton method. In the following lines, we demonstrate that the derivative f 0 is of a constant sign, which implies a unique solution for (15), if it exists.

The derivative of the function f with respect to λ is as follows f 0 λ = r 0 λ1 , p rλ + 1 2 ,1

It is straightforward to prove that r 0 λ = λ , p rλ so we have

f 0 λ = λ , p1 , p rλrλ + 1 2 ,1 (16)
Let us call P 0 is the projection of P on (CB). We have then λ , p = , HP 0 l 1 , p = , BP 0 l To prove the constant sign of f 0 , we just have to compare the two distance products CB:HP 0 :BP 0 (the numerator of the In our implementation, we use a Newton method to find λ 0 solution of (15). Starting from a central position (λ = 1), after 3 or 4 iterations, the difference between two successive estimations of λ do not exceed 10 ,6 . The convergence is extremely fast.

B. PROLATE SPHEROIDAL COORDINATES

We define a center O and two "focal points" F 1 and F 2 , F 2 being at the same distance δ from O as F 1 , but in the opposite direction. A prolate sphere is defined to have a constant radius λ (dimensionless number), a point in this prolate sphere is defined fixing two angles: elevation θ and azimuth φ. From these three parameters, the cartesian coordinates (x, y, z) of this pont are calculated using the following formulae:

x = δ sinhλ sinθ cosφ y = δ sinhλ sinθ sinφ ( 17)

z = δ coshλ cosθ
Conversely, knowing the cartesian coordinates (x, y, z), it is possible to compute the prolate spheroidal parameters (λ,θ,φ) using the equations:

r 1 = q x 2 + y 2 + z ,δ 2 r 2 = q
x 2 + y 2 + z + δ 2 λ = acosh 0 @ r 1 + r 2 2:δ 1 A θ = acos 0 @ r 1 ,r 2 2:δ 1 A (18) φ = atan2 0 @ x δ sinhλ sinθ ; y δ sinhλ sinθ 1 A Analogously to our system, F 1 would be C, O would be B, F 2 would be C 0 (figure 7) and λ would be R. In the end, the transformation from (R,θ,φ) coordinates to (X,Y,R) would be expressed as in (2).

Figure 2 .

 2 Figure2. The top three frames illustrate the orthogonal motions described in the text. Bottom, point P is transformed in P 1 by the first motion (radial motion), P 1 is transformed in P 2 by the second (apico-basal rotation) and P 2 is transformed in P 3 by the third one (elevation).

  Figure 3. P 1 belongs to a lateral wall, P 2 is close to the apex, both points are animated with motion 1. Solid arrows describe the motion in cylindrical coordinates (contraction towards the apico-basal axis) and dotted arrows describe the motion in spherical coordinates (contraction towards the center C). Motion 1 effectively describes a contraction for P 1 , but not for P 2 .

Figure 4 .

 4 Figure 4. 3D-planispheric geometry is a combination of both cylindrical and spherical geometries.

FigureFigure 5 .

 5 Figure5. the cartesian coordinates (x,y,z) of point P are converted into polar coordinates in the 3D-planispheric image, the depth R is the distance from the point H P in the cartesian image, the position (X,Y) in the plane is defined with the two angles θ and φ, like in the 2D-planispheric mapping.

Figure 7 .

 7 Figure7. A representation of ten surfaces of constant R, on the left, in the 3PC system, on the right, in the prolate spheroidal. The dashed lines show different curves of constant θ: on one, a point P and its associated center H P . Figure from(Declerck et al,1997).

  Figure8. The matching criterion measures the sum of distances between f (P l , t n ), the estimated point P l at time t n and Q n;l , the closest feature point of Image n, for all existing P l and t n .

Figure 9 .

 9 Figure 9. The template and three points defining the coordinate system. By the transformation S, they are deformed to match the shape of the left ventricle of the patient.

  anterior vector : w T = u T v T For the patient's case -center : C P = SC T -base : B P = SB T -lateral point : L P = SL T -apicobasal vector : u P = C P B P kC P B P k -septo-lateral vector :v P =

Figure 10 .

 10 Figure 10. The surface S is represented in 3PC. In this geometry, the location of the closest point Q to P of surface S depends on σ r . Q belongs to the curved segment [Q 0 Q ∞ ], where Q 0 and Q ∞ are the closest point to P for σ r = 0 and ∞ respectively. Figure reprinted with kind permission of IEEE Trans. on Med. Imag.

Figure 11 .

 11 Figure 11. Intensity profiles along 2 different radii starting from the center of the image. The black dots show on both curves the edges detected with the Canny-Deriche filter. On curve 2, the noncardiac edge is rejected. Figure reprinted with kind permission of IEEE Trans. Med. Imag.

Figure 12 .

 12 Figure 12. Edges (in white) automatically extracted and filtered from images at four different times (one time per row). On each row, we see central slices resampled from the 3D image by rotation around the apex-base axis, as shown on drawing on the left: this makes easier the display of the myocardial structure.

Figure 13 .

 13 Figure13. View of the LV from the apex. The trajectories of some points are drawn over the cycle: they are smooth and periodic (see the zoomed area on the right). Figure from(Declerck et al,1997).

  Figure13. View of the LV from the apex. The trajectories of some points are drawn over the cycle: they are smooth and periodic (see the zoomed area on the right). Figure from(Declerck et al,1997).

Figure 16 .Figure 17 .

 1617 Figure 16. The surface of the heart (the septum is on the left) colorised with the parameters values at end systole, for a normal case (VANN, on the left) and for a pathological case (BOJO, on the right): top, the radial motion, middle, the elevation and bottom, the apicobasal rotation. The colorscales under the images indicate increasing values from left to right. Figure from(Declerck et al,1997).

  BC; HP = θ CH = 1 ,cosθ:CB Let us define λ as follows: CH = λ:CB = λ:l Because λ is supposed to be equal to 1 ,cosθ, λ 2 0= u:HC + u:CP = l: ,λ+ x:u x + y:u y + z:u z On the other hand, :u x 2 + y ,λ:u y 2 + z ,λ:u z 2 and p = u:CP we thus have cosθ = , u:HP r and the constraint (1) can be rewritten, after some calculation, f a function of λ f λ = rλ + p rλ + 1 ,λ

Figure 18 .

 18 Figure 18. Two possible configurations for P 0 .

, and PH P is nearly orthogonal to CB: around the base, the

3PCsystem is close to the cylindrical one (figure7).

ACKNOWLEDGEMENTS

We want to thank Pr. Michael L. Goris for his precious advices we had and the team of the Stanford University Hospital (California, USA) for providing us with the images.

We give also special thanks to Dr. Eric Bardinet and Pr. Mike Brady for the constructive discussions and comments about the project.

This work was partially supported by regional grant of the Région Provence Alpes Côte d'Azur (doctoral research contract).