Eric Bardinet

Laurent D Cohen

Nicholas Ayache

Tracking and motion analysis of the left ventricle with deformable superquadrics

Keywords: Parametric models, tracking, left ventricle motion, quantitative analysis of motion

We present a new approach to analyse the deformation of the left ventricle of the heart based on a parametric model that gives a compact representation of a set of points in a 3-D image. We present a strategy for tracking surfaces in a sequence of 3-D cardiac images. Following tracking, we then infer quantitative parameters which characterize: left ventricle motion, volume of left ventricle, ejection fraction, amplitude and twist component of cardiac motion. We explain the computation of these parameters using our model. Experimental results are shown in time sequences of two modalities of medical images, nuclear medicine and X-ray computed tomography (CT). Video sequences presenting these results are on the CD-ROM.

INTRODUCTION

The analysis of cardiac deformations has given rise to a large amount of research in medical image understanding. Indeed, cardiovascular diseases are the primary cause of mortality in developed countries. Various imaging techniques [START_REF] Acharya | Biomedical imaging modalities: a tutorial[END_REF] allow the acquisition of dynamic sequences of 3-D images (3-D+T) during a complete cardiac cycle (contraction and dilation). These images are well-suited to studying the behaviour of the cardiac system since they visualize how the heart wall deforms. Processing these images opens numerous fields of applications, such as detection and analysis of pathologies.

Advanced techniques of 3-D imagery, such as nuclear medicine and X-ray computed tomography (CT) provide ever increasing resolution in space and time. Consequently, the data available to the radiologist are becoming larger. However, to establish a reliable and fast diagnosis, the physician needs models that are defined by only a small number of characteristic quantities. A parametric deformable model allows the representation of a dynamic set of points by a

• The left ventricle reconstruction was done with generic deformable surface models [START_REF] Ayache | Building highly structured volume representations in 3D medical images[END_REF][START_REF] Cohen | On active contour models and balloons[END_REF][START_REF] Leitner | EMBS '91[END_REF]Cohen et al., 1992b;[START_REF] Cohen | Finite element methods for active contour models and balloons for 2-D and 3-D images[END_REF][START_REF] Ayache | Computer-Integrated[END_REF] but also with surface models dedicated to the left ventricle shape (Duncan et al., 1991b;[START_REF] Amini | Bending and stretching models for LV wall motion analysis from curves and surfaces[END_REF][START_REF] Clarysse | 3-D boundary extraction of the left ventricle by a deformable model with a priori information[END_REF]. • The left ventricle tracking was also studied with generic deformable surface models [START_REF] Pentland | Recovery of nonrigid motion and structure[END_REF][START_REF] Ayache | Active Vision[END_REF][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] and with the help of curvature information [START_REF] Amini | Nonrigid motion models for tracking the left ventricular wall[END_REF]Duncan et al., 1991a;Cohen et al., 1992a;[START_REF] Friboulet | Assessment and visualization of the curvature of the left ventricle from 3-D medical images[END_REF][START_REF] Benayoun | Adaptive meshes and nonrigid motion computation[END_REF][START_REF] Benayoun | Dense non-rigid motion estimation in sequences of 3D images using differential constraints[END_REF]. Four dimensional models have been proposed by [START_REF] Shi | Shape-based 4D left ventricular myocardial function analysis[END_REF] and [START_REF] Nastar | Frequency-based nonrigid motion analysis: application to 4 dimensional medical images[END_REF] and the exploitation of temporal constraints was studied by [START_REF] Meyer | Tracking myocardial deformation using spatially-constrained velocities[END_REF] and [START_REF] Mceachen | A recursive filter for phase velocity assisted shapebased tracking of cardiac non-rigid motion[END_REF]. • Finally, the extraction of parameters which capture the overall deformation was presented by Duncan et al. (1991a), [START_REF] Shi | A model-based integrated approach to track myocardial deformation using displacement and velocity constraints[END_REF], [START_REF] Benayoun | Dense non-rigid motion estimation in sequences of 3D images using differential constraints[END_REF] and [START_REF] Nastar | Frequency-based nonrigid motion analysis: application to 4 dimensional medical images[END_REF]. • In some images, some sparse 'anchor points' can be produced within the image to help the tracking process. This the case of MRI-SPAMM images, and a number of studies take advantage of their properties [START_REF] Amini | MR physics-based snake tracking and dense deformations from tagged cardiac images[END_REF][START_REF] Guttman | Tag and Contour Detection in Tagged MR Images of the Left Ventricle[END_REF][START_REF] Young | Deformable Models for Tagged MR Images: Reconstruction of 2-and 3-Dimensional Heart Wall Motion[END_REF][START_REF] Kumar | Automatic Tracking of SPAMM Grid and the Estimation of Deformation Parameters from Cardiac MR Images[END_REF][START_REF] Park | Analysis of left ventricular wall motion based on volumetric deformable models and MRI-SPAMM[END_REF].

In a previous article (Bardinet et al., 1996a), we introduced a parametric deformable model based on a superquadric fit followed by a free form deformation (FFD). The advantage of parametric deformable models like superquadrics is the small number of parameters needed to describe a shape combined with a better robustness in the presence of noise or sparse data. Also, at the expense of a reasonable number of additional parameters, FFDs provide a close fit and a volumetric deformation estimation.

In the present article, we first give a summary of our segmentation algorithm, specific to cardiac images, in Section 2. In Section 3, we briefly summarize the results of Bardinet et al. (1996a) on the parametric model, necessary for a good understanding of this article. In Section 4, we present an approach to track surfaces with this model in a sequence of 3-D images, and give experimental results for tracking the deformation of the left ventricle in two different kinds of 3-D medical images. In Section 5, we explain how to infer, from the parametric reconstruction, a number of quantitative parameters useful to characterize the left ventricle motion. We demonstrate the feasibility of the approach on two kinds of temporal sequences of 3-D images. We believe that the results are sufficiently promising to initiate a thorough clinical validation.

SEGMENTATION OF CARDIAC IMAGES

We studied two different kinds of images.

• Nuclear medicine data, the SPECT sequence, with eight successive time frames during one cardiac cycle. Each image is a volume of 64 × 64 × 64 voxels (volume of the voxel: 1.446 mm 3 ), describing a human heart. • X-ray CT data, the DSR (dynamic spatial reconstructor) sequence, with 18 successive time frames during one cycle. Each image is a volume of 98 × 100 × 110 voxels (volume of the voxel: 0.926 mm 3 ), describing a canine heart.

The original 3-D images are visualized as a series of 2-D cross-sections in Figures 1 and3. We first have to extract a set of points belonging to the endocardium (the inner surface of the left ventricle) and/or to the epicardium (the outer surface of the left ventricle). These points will be then approximated by a deformable superellipsoid in the next section. In DSR images, a single threshold is sufficient to isolate the left ventricle cavity, and therefore extract the endocardium. On the other hand, the epicardium is not sufficiently well-contrasted to be robustly extracted. In the SPECT images, it is possible to extract grossly the epicardium and endocardium surfaces with a threshold based on the histogram of the intensities (Goris and Bertille, 1992). Then, with the help of mathematical morphology operators [START_REF] Serra | Image Analysis And Mathematical Morphology[END_REF][START_REF] Höhne | Interactive 3D segmentation of MRI and CT volumes using morphological operations[END_REF], we automatically smooth and isolate both surfaces.

Therefore, in the following, we assume that we have extracted the points belonging respectively either to the endocardium (DSR and SPECT images) and/or to the epicardium (SPECT images only). Although extracting surfaces from the SPECT sequence seems to be a difficult task, let us point out that we ran experiments successfully on six SPECT sequences with the same segmentation process. Results of the segmentation processes on the two data sequences are presented in Figures 2 and4 for data at time step 1, and in Figures 6 and8 for one cross-section over time during the sequence. For a correct estimation of the quality of those segmentations, we superimposed the segmented surface(s) on the image.

A PARAMETRIC MODEL TO FIT 3-D DATA

In this section, we briefly describe the deformable model that we use to represent the inner and outer surfaces of the left ventricle [a detailed presentation is in Bardinet et al. (1996a)].

We refine our superquadric model using a parametric deformation. More precisely, for a given set of 3-D points (we have seen in the previous section how get from the 3-D images of the heart a number of 3-D points belonging to the inner or outer surfaces of the left ventricle), we first fit 3-D data with a superellipsoid, and then refine this crude approximation using free form deformations (FFDs).

Fitting 3-D data with superquadrics

Superquadric shapes have been widely used in vision and graphics. In computer vision, their first use is due to [START_REF] Pentland | Recognition by parts[END_REF], followed by [START_REF] Solina | Recovery of parametric models from range images: the case for superquadrics with global deformations[END_REF] who used superellipsoids to approximate 3-D objects. The goal of the algorithm is to find a set of parameters such that the superellipsoid best fits the set of data points. Superquadrics form a family of implicit surfaces obtained by extension of conventional quadrics. Superellipsoids centered at the origin and with principal axes corresponding to the reference axes are defined by the implicit equation:

  x a 1 2/ǫ 2 + y a 2 2/ǫ 2 ǫ 2 /ǫ 1 + z a 3 2/ǫ 1   ǫ 1 /2 = 1, (1)
which involves five independent parameters. To generate a superellipsoid centered at an arbitrary location and with an arbitrary orientation of its principal axes, we must add the six parameters of a rigid displacement, which makes 11 parameters required for an arbitrary superellipsoid. Suppose that the data we want to fit with the superellipsoid are a set of 3-D points (x i , y i , z i ), i = 1, • • • , N . Since a point on the surface of the superellipsoid satisfies F = 1, where F is the function defined by Equation (1), we find the minimum of the following energy:

E(A) = N i=1 [1 -F(x i , y i , z i , A)] 2 , ( 2
)
where A denotes the set of 11 parameters needed to describe the superellipsoid [see Bardinet et al. (1996a) for details about the minimization technique used and also for a geometric interpretation of this energy].

Refinement with free form deformations

We now refine this parametric representation of the 3-D data, using a global volumetric deformation called FFD. This is a tool devoted to the deformation of solid geometric models in a free-form manner [see [START_REF] Sederberg | Free-form deformation of solid geometric models[END_REF]]. The main interest of FFDs is that the resulting deformation of the object, although potentially complex, is defined by a small number of points. This characteristic feature allows us to represent voluminous 3-D data by a model defined by a relatively small number of parameters (11 parameters + 125 points).

Definition of FFDs

FFDs were introduced by [START_REF] Sederberg | Free-form deformation of solid geometric models[END_REF] in computer graphics and have been used to solve matching problems by [START_REF] Szeliski | Matching 3-d anatomical surfaces with non-rigid deformations using octree-splines[END_REF]. An FFD is a mapping from R 3 to R 3 , defined by the tensor product of trivariate Bernstein polynomials. The principle of FFDs is as follows: the object to be deformed is embedded in a 3-D box.

Inside this box, a volumetric grid of points is defined, which links the box to the object (by the trivariate polynomial which defines the deformation function):

X = l i=0 m j=0 n k=0 C i l C j m C k n (1 -s) l-i s i ×(1 -t) m-j t j (1 -u) n-k u k P i jk , (3)
where P i jk denotes the volumetric grid of control points, (s, t, u) denotes the local coordinates of the object points in a frame defined by the box of control points and (l, m, n) denotes the degrees of the Bernstein polynomials. This can be written in a matrix form: X = BP, where B is the deformation matrix N d × N p (N d is the number of points on the discretized superellipsoid and N p is the number of control points of the grid), P is a matrix N p × 3 which contains the coordinates of the control points and X is a matrix N d × 3 with the coordinates of the model points. The box is then deformed by the displacement of its lattice, and the position of a point of 2 × 2 × 2 = 8 243 3 × 3 × 3 = 27 244 4 × 4 × 4 = 64 265 5 × 5 × 5 = 125 429 6 × 6 × 6 = 216 751 the deformed object is computed [see Bardinet et al. (1996a) for details].

The inverse problem

We need to solve the inverse problem: first compute a initial displacement field δX between the superellipsoid and the data and then, after having put the superellipsoid in a 3-D box, search for the deformation δP of this box which will best minimize the displacement field δX:

min δP BδP -δX 2 (4)
This is illustrated in Figure 5.

In practice, we use an iterative two-step algorithm [see Bardinet et al. (1996a) for details], similar to the formulation of B-spline snakes with auxiliary variables, as described in [START_REF] Cohen | Auxiliary variables and two-step iterative algorithms in computer vision problems[END_REF].

The current displacement field δX is computed after a set of matches is established between the image data points and the current position of the model (which is a deformed superellipsoid after the first iteration). In our implementation, a data point is matched with the closest model point, but other strategies, involving for instance the local curvature, could be successfully applied (Duncan et al., 1991a;[START_REF] Amini | Nonrigid motion models for tracking the left ventricular wall[END_REF]Cohen et al., 1992a;[START_REF] Ayache | Active Vision[END_REF][START_REF] Shi | Shape-based 4D left ventricular myocardial function analysis[END_REF][START_REF] Benayoun | Adaptive meshes and nonrigid motion computation[END_REF]Feldmar and Ayache, 1994;[START_REF] Benayoun | Dense non-rigid motion estimation in sequences of 3D images using differential constraints[END_REF].

We use boxes of size 5 × 5 × 5 for data sets of ∼6,000 points, and the number of iterations is between 10 and 30. Therefore we get an average compression ratio of 46 [computed as (6000 × 3)/(125 × 3 + 11), where 11 is the number of parameters needed for the superellipsoid (compare with Subsection 3.1)]. This number of parameters provides a reasonable trade-off between approximation quality and computation time. Table 1 presents typical computation times of the FFD with 20 iterations and different sizes for the box of control points (on a DEC Alpha 300).

Simultaneous deformation of two surfaces

An essential feature of this algorithm is that FFD is a volumetric deformation. This means that several objects can 4) is done simultaneously based on the union of the two displacement fields.

Figure 5 shows the result of the algorithm for the reconstruction of the epicardium and the endocardium from the SPECT data with two superellipsoids deformed first by two independent FFDs, and then by a single common FFD (in this last case, the unique box of control points is the largest one, i.e. the one computed for the epicardium). The corresponding approximation errors for the SPECT data, computed as the averages of errors between the data and the final model (projecting each 3-D data point onto the deformed model), are presented in Table 2. One can see that using two FFDs for the two surfaces leads to a better quality of approximation. Left column, each model is computed independently; middle column, the two models are computed with one FFD; right column, computation with only one FFD leads to a slight loss of the approximation precision.

On the other hand, using only one FFD allows to divide the number of parameters by two, yielding a larger compression of the information needed for the description of the parametric model.

TRACKING OF THE LEFT VENTRICLE

Having extracted 3-D points corresponding to the wall of the left ventricle in a time sequence of images, we show in this section how we use our parametric model to make an efficient tracking of the LV wall in a sequence of 3-D images.

Tracking strategy: recursive representation with a unique deformation

Our strategy is the following: the complete model (superellipsoid + FFD) is fitted only to the surfaces of the left ventricle extracted in the first image. Then, a unique FFD (which is now a deformed box of control points) is used with a new displacement field computed between the previous model and the data of the next image, and so on for the complete sequence. Finally, the models at any time are defined by only one FFD applied to a single original superellipsoid.

Following the formulation of Subsection 3.2.2, the model at time n, X n , can be written:

X n = B 0 (P 0 + δP * 0 + • • • + δP * n-1 ) , (5)
where δP * i is the increment in the position of control points at iteration i which minimizes the error defined in Subsection 3.2.2 between the deformed model and the extracted surface(s) of the left ventricle.

One can notice that this strategy requires the inversion of a linear system of the form X = BP only once [by computing the pseudoinverse of the matrix with a singular value decomposition, see Bardinet et al. (1996a) for details], thus tracking of a sequence is quite efficient.

Experimental results on cardiac images

We present in this section results of the tracking algorithm on the cardiac images described in Section 2.

Figure 6 shows the dynamic SPECT sequence on a particular cross-section (z = 32), with the segmented and reconstructed surfaces superimposed (see Section 2). It also shows the reconstruction of these two surfaces using either one or two models as explained in Subsection 3.2.3. Figure 7 presents the same results for the 3-D rendering of the surface.

Figure 8 shows the dynamic DSR sequence on a crosssection, first the original image, then the segmented endocardium surface superimposed (see Section 2) and also the reconstruction of this surface with the parametric model. Figure 9 depicts the same result for the 3-D rendering of the surface for three time steps of the sequence. Table 3 shows the CPU time for each complete time sequence (on a DEC Alpha 300).

All these results are also presented in the video attached to this article (see Appendix C).

QUANTITATIVE ANALYSIS OF LEFT VENTRICLE DEFORMATION

The reconstruction and representation of a time sequence of surfaces by a sequence of parametric models allows to visualize the estimation of the deformation through time.

More precisely, the parametric representation provides a way to determine the motion field on the cardiac wall. This motion field can then be used to extract some characteristic parameters useful for the diagnosis.

Let us here point out that these parameters are deduced from a geometrical representation of the left ventricle motion. Thus, the question we have to keep in mind is the following: what kind of parameters are of interest for the physician? There are parameters which he already uses for diagnosis, and he is familiar with, so that he can assess the practicability of the model. In the area of cardio-vascular diseases, and especially in the study of the cardiac muscle, an important parameter is the volume of the left ventricle, as well as its temporal evolution, which permits to compute the ejection fraction (ratio characterizing the 'pumping' nature of the heart). This is the content of Subsection 5.1.

We also present other types of parameters. These ones are not yet used for diagnosis, and a complete clinical validation is clearly required to assess their usefulness.

In Subsection 5.2, we present parameters that can be deduced from the trajectories of parametric points of the model during a cardiac cycle: amplitude and twist component of the motion, temporal evolution of the left ventricle wall thickness.

Volume evolution

To evaluate the ejection fraction, we need to compute the temporal evolution of the left ventricle volume. We first calculate an explicit form of the volume of a superellipsoid, depending on its parameters, and then show how to calculate the volume of a polyhedric region bounded by a mesh of vertices, using the discrete form of the Gauss integral theorem.

Volume of a superellipsoid

Let S be a superellipsoid surface defined by the following implicit equation:

x a 1 2 ǫ 2 + y a 2 2 ǫ 2 ǫ 2 ǫ 1 + z a 3 2 ǫ 1 = 1. ( 6
)
Its volume can be explicitly computed by using the following formula:

V = 2 a 1 a 2 a 3 ǫ 1 ǫ 2 β ǫ 2 2 , ǫ 2 2 + 1 β ǫ 1 2 , ǫ 1 + 1
(see Appendix A for details).

Volume inside a mesh

Considering an oriented mesh defined by a set of points and a set of faces, the volume inside the mesh can be written:

V = K i=1 V i = 1 3 K i=1 (O M, N ) S i , ( 7
)
where K is the number of faces of the mesh, O is a reference point, M is a point on the mesh and N is the normal vector at M, pointing outwards. This expression can be rewritten:

V = 1 6 K i=1 (OG, N d ) ,
where N d and G denote respectively the normal vector and the barycenter of a face (see Appendix B for details).

Discussion

We applied this calculation of epicardium and endocardium volumes to the sequences of both data points and parametric models obtained in the previous sections. Once we have computed the values of the volume along a cardiac cycle, we can easily obtain the ejection fraction [calculated precisely as:

Vd-Vc Vd , with Vd volume at dilation (the end of diastole), V c volume at contraction (the end of systole), see for example [START_REF] Davis | Assessment of left ventricular ejection fraction from technetium-99m-methoxy isobutyl isonitrile multiple-gated radionuclide angiocardiography[END_REF]]. The results presented in Figure 11 show the following.

• The evolution of the volume has the expected typical shape found in the medical literature [START_REF] Goris | A Colour Atlas of Nuclear Cardiology[END_REF] and shown in Figure 10. Moreover, the estimation of the ejection fraction on our example gives a value of 68%, which is in the range of expected values [START_REF] Goris | A Colour Atlas of Nuclear Cardiology[END_REF][START_REF] Davis | Assessment of left ventricular ejection fraction from technetium-99m-methoxy isobutyl isonitrile multiple-gated radionuclide angiocardiography[END_REF]. • The volumes computed on the data points or on the deformable superellipsoids are almost equal. The relative average absolute error along the cycle is 0.42%. This proves that our model is robust with respect to the volume estimation. Of course, the ejection fraction is also obtained with a very small relative error (0.19%). • The volume evolution found for initial superellipsoid models before FFD, has also a very similar shape, as seen in Figure 11. However, there is a size ratio due to the over-estimation of the volume before the FFD. This ratio is almost constant through time, which makes it possible to get a good estimate of the ejection fraction directly from the initial model. This proves that the superellipsoid model provides a good global estimate of the shape. Also, the volume of the superellipsoid can be obtained analytically from its set of parameters without the previous discrete approximation (see Subsection 5.1.1).

Trajectories

Listing the successive positions of a parametric point of the deformed surface model along the time sequence, we obtain the trajectory followed by this point.

Decomposition into cylindrical coordinates

The decomposition of the model points trajectories into cylindrical coordinates (see Figure 12 and temporal evolution of the wall thickness.

  x y z   →     ρ = x 2 + y 2 θ = arcos x x 2 + y 2 z = z     (8)
We will consider, as frame of reference of the cylindrical representation, the frame of reference of the superellipsoid fitted to the segmented data at the first step of the time sequence.

To measure the twist, we compute the difference of the θ parameters for the two points that represent the same parametric point during the contraction.

The wall thickness is computed as the difference of the ρ parameters for two corresponding parametric points on the epicardium and the endocardium [correspondences can be established by different criteria: in a preliminary study, we used the parametric values on the original superellipsoids, before the application of FFDs, see [START_REF] Bardinet | Analyzing the deformation of the left ventricle of the heart with a parametric deformable model[END_REF] for details, but this will be improved in a future work].

Synthetic example

To get a precise idea of the capability of the model in recovering a twist motion, we studied a synthetic example.

First, we define a sequence of parametric surfaces representing a synthetic heart during a complete cardiac cycle (17 successive time frames). All the models of this sequence are obtained with the following implicit equation: where a 1 = 10, a 2 = 12, a 3 = 14 and c 1 and c 2 are periodic temporal functions to simulate a gross beating motion (a more sophisticated model would have included a translation of the origin and a rotation of the principal axes). The complete sequence is computed with the following 17 values of c 1 : 3.3, 2.9, 2.5, 2.3, 2.1, 1.9, 1.7, 1.5, 1.7, 1.9, 2.1, 2.3, 2.5, 2.7, 2.9, 3.1, 3.3 (same values for c 2 ).

x a 1 2 + y a 2 2 + z a 3 2 3 -c 1 x a 1 2 + c 2 y a 2 2 z a 3 3 = 1,
Starting from this sequence, we add a global axial twist around the z-axis of the models, as defined by [START_REF] Barr | Global and local deformations of solid primitives[END_REF]:

θ = f (z) C θ = cos(θ), S θ = sin(θ)    X = xC θ -y S θ Y = x S θ + yC θ Z = z
f is defined as follows: the angle θ varies linearly from α on the apex to -α on the base, with α varying along the sequence from 0 o up to 30 o , then back to 0 o . A very recent work (Wakes et al., 1996) presents a cardiac motion simulator, which could yield a more sophisticated synthetic example.

We then try to recover the known trajectories using our tracking strategy. Recovering a twist motion is a difficult problem, since the displacement of each point has a tangential component which cannot be easily detected. Figure 13 presents the recovered trajectories compared with the original ones.

One can see (and this is particularly visible in the attached video) that the tracking method provides results which are close to the original trajectories, and catch the twist component of the motion. This is probably due to the discrete nature of the synthetic example, which forces the 'closest' match on the second surface to include a tangential component. With real data, either tags or curvature attributes would be required to obtain a similar result.

Figure 14 represents the values of this twist, where a colour has been associated to the twist value for each model point, on the original synthetic data and for the models recovered by the tracking method.

Real images

Figure 15 show the estimated trajectories on the epicardium surface during the systolic stage for the SPECT data. These are visualized from different viewpoints to appreciate better the motion, and also in the video attached to this article (see Appendix C).

The pointwise tracking of the deformation allows an evaluation of the motion field during the sequence. The visualization of these displacements using different colour values on the surface clearly shows those areas on the ventricle where the displacements are small (see Figure 16).

CONCLUSION

We have presented a new approach to analyse the deformation of the heart left ventricle with a parametric model. It is based on a parametric model that gives a compact representation of a set of points in a 3-D image. We have presented a method which efficiently tracks the left ventricle wall in a sequence of 3-D images during a cardiac cycle using this model and studied synthetic and real data. The model is able to track simultaneously the endocardium and the epicardium, since it is a volumetric deformation.

It is then possible to compute the variation of the volume during a cardiac cycle and the ejection fraction. We showed on a synthetic example that it might also be possible to estimate the twist component in the deformation of the ventricle, provided that some local features allow for the recovery of the tangential component of the motion.

Future directions will include the clinical evaluation of this work, and also its adaptation to MRI-SPAMM images where we could use our model to track and interpolate a smooth deformation between tagged points.

APPENDICES A. CALCULATION OF THE VOLUME OF A SUPERELLIPSOID

Let S be a superellipsoid surface defined by the following implicit equation:

x a 1 2 ǫ 2 + y a 2 2 ǫ 2 ǫ 2 ǫ 1 + z a 3 2 ǫ 1 = 1. (9)
An explicit parameterization of S is given by:

   x = a 1 cos ǫ 1 η cos ǫ 2 ω y = a 2 cos ǫ 1 η sin ǫ 2 ω z = a 3 sin ǫ 1 η , -π 2 ≤ η ≤ π 2 -π ≤ ω < π
Due to the symmetry of the surface in relation to the 3 axes of the coordinate system, the computation of the volume inside S can be made as follows:

V = 2 a 3 0 A(z) dz, ( 10
)
where A(z) is the area of a slice among the z-axis. For η = 0, the implicit equation of the corresponding slice is:

x a 1 2 ǫ 2 + y a 2 2 ǫ 2 = 1
This leads to:

A(0) = 4 a 2 a 1 0 y(x) d x = 4 a 2 a 1 0 1 - x a 1 2 ǫ 2 ǫ 2 2 d x Setting: X = x a 2 ǫ 2
2 , it becomes:

A(0) = 2 a 1 a 2 ǫ 2 β ǫ 2 2 , ǫ 2 2 + 1 ,
where β(x, y) denotes the Beta function (Euler's integral of the first kind):

β(x, y) = 1 0 t x-1 (1 -t) y-1 dt
Finally, to calculate the volume inside S using Equation 10, we have to write a 1 and a 2 as functions of z. From the implicit definition of S (Equation 9), we deduce:

a 1 (z) = a 1 1 - z a 3 2 ǫ 1 ǫ 1 2 , a 2 (z) = a 2 1 - z a 3 2 ǫ 1 ǫ 1 2 Therefore: V = 2 a 3 0 2 ǫ 2 β ǫ 2 2 , ǫ 2 2 + 1 a 1 (z) a 2 (z) dz Setting: Z = z a 3 2 ǫ 1
, it becomes:

V = 2 a 1 a 2 a 3 ǫ 1 ǫ 2 β ǫ 2 2 , ǫ 2 2 + 1 β ǫ 1 2 , ǫ 1 + 1
Note that for a sphere (ǫ 1 = ǫ 2 = 1, a 1 = a 2 = a 3 = R), the previous formula gives:

V = 2 R 3 β 1 2 , 3 2 β 1 2 , 2 = 4 3 π R 3

B. CALCULATION OF THE VOLUME INSIDE A MESH

Let D be a region of space with bounding surface ∂ D; the unit normal vector n to ∂ D is drawn outwards. Then: known as the Gauss integral theorem, where div(.) denotes the divergence of a vector field and (., .) denotes the scalar product in R 3 .

Let O be a reference point. Using this theorem, we can write:

D div(O M) d x dy dz = ∂ D (O M, n) d S,
where M is a point on ∂ D. Now, it is obvious that the value of the divergence of the vector field O M is 3. Therefore, the Gauss integral thereom yields an expression of the volume inside a domain D:

V = D d V = 1 3 ∂ D (O M, n) d S ( 12
)
Considering an oriented mesh defined by a set of points and a set of faces, the volume inside the mesh can be written, using the previous formula:

V = K i=1 V i = 1 3 K i=1 (O M, N ) S i , ( 13
)
where K is the number of faces of the mesh. We now have to define the normal vector of a face. Assuming that the face is defined by three points A 1 , A 2 and A 3 , the normal vector can be written:

N d = (A 1 -A 2 )∧(A 1 -A 3 ) = A 1 ∧ A 2 + A 2 ∧ A 3 + A 3 ∧ A 1 ,
where (. ∧ .) denotes the vector product. Note that the norm of this normal vector is linked to the surface S of the face as follows: N d = 2S . One can easily generalize this formula for a face defined by a set of points A i , i = 1..M:

N d = M-1 i=1 (A i ∧ A i+1 ) + A M ∧ A 1
This expression can be seen as the average of the normal vectors of the decomposition of the face into triangular subfaces.

Let suppose that the mesh is composed of triangular faces defined by A 1 , A 2 and A 3 . The volume V i of the tetrahedron O A 1 A 2 A 3 is:

V i = S.H 3 ,
where H is the height. Let O be the origin and G be the barycenter of a face. Then:

H = (OG, N d ) N d
And therefore: V i = (OG, N d ) 6 This last expression is still true for a face defined by M points. Finally, Equation 13 can be rewritten:

V = 1 6 K i=1 (OG, N d )

C. CONTENT OF THE VIDEO

The video is consisting in six sequences:

1. Nuclear medicine data (SPECT image). The first image of the sequence is visualized as a series of 2-D crosssections (transverse slices). See Section 2 for details (0-20 s).

2. Segmentation of the SPECT sequence. See Section 2 for details (20-40 s).

3. Tracking of the epicardium (mesh) and the endocardium (rendered surface) in the SPECT sequence; on the left, the segmented surfaces, on the right, the reconstructed models. See Sections 2 and 4.2 for details (40 s-1 min).

4. On the left, trajectories of the node points, on the right, velocity field in the SPECT sequence. See Subsection 5.2.3 for details (1-1 min 20 s).

5. Recovery of a twist motion on a synthetic example (1). Data with trajectories. See Subsection 5.2.2 for details (1 min 20 s-1 min 40 s).

6. Recovery of a twist motion on a synthetic example (2).

Trajectories on the data and on the models computed with tracking method. See Subsection 5.2.2 for details (1 min 40 s-2 min).

Figure 1 .

 1 Figure 1. 3-D image of the left ventricle -SPECT image (the order of sections reads from left to right and from top to bottom). Dynamic presentation in the video.

Figure 2 .

 2 Figure 2. Segmentation of the epicardium and the endocardium (external and internal surfaces of the left ventricle) on the SPECT image (the order of sections reads from left to right and from top to bottom). Dynamic presentation in the video.

Figure 3 .

 3 Figure 3. 3-D image of the left ventricle -DSR image (the order of sections reads from left to right and from top to bottom).

Figure 4 .

 4 Figure 4. Segmentation of the endocardium (internal surface of the left ventricle) on the DSR image (the order of sections reads from left to right and from top to bottom).

Figure 5 .

 5 Figure 5. Deformation of two surfaces. Top left: data: epicardium (mesh) and endocardium (rendered surface) described respectively with 4500 and 1500 points. Top right: two fitted superellipsoids and the two initial boxes of control points. Bottom left: final models obtained with two FFDs (each defined by 11 parameters + 125 points). Compression ratio: 23. Bottom right: simultaneous deformation of the two superellipsoids with only one FFD (defined by 11 parameters + 125 points). Compression ratio: 46.

Figure 6 .

 6 Figure 6. Segmentation and representation of the epicardium and the endocardium for a time sequence (3-D+T) of the left ventricle during the cardiac cycle -SPECT image; visualization of a cross-section over time. Top: segmented surfaces superimposed on the image. Middle: approximation by 2 independent deformable superellipsoids (2 FFDs). Bottom: approximation by 2 coupled deformable superellipsoids (1 FFD).

Figure 7 .

 7 Figure 7. Time sequence of the epicardium (mesh) and the endocardium (rendered surface) (t = 1, 3, 5) SPECT images. Left: isosurfaces obtained by data segmentation (4500 + 1500 points). Right: representation by two parametric models (2× (11 parameters + 125 points)). Dynamic presentation in the video.

Figure 9 .

 9 Figure 9. Time sequence of the endocardium (t = 1, 8, 13) DSR images. Left: isosurfaces obtained by segmentation (10000 points). Right: reconstruction by the parametric model defined by 11 parameters + 125 points (DSR). Compression ratio: 77.

Figure 10 .

 10 Figure 10. Typical shape of the temporal evolution of ventricular volume. The curve is divided into 3 parts: 1. Fast contraction 2. Fast dilation 3. Slow dilation.

Figure 11 .

 11 Figure 11. Cardiac volumes during the cardiac cycle. Left: volumes of the epicardium (SPECT). Middle: volumes of the endocardium (SPECT). Right: volumes of the endocardium (DSR). Top row: volumes computed on the data points. Middle row: volumes of the superellipsoids. Bottom row: volumes of the deformable superellipsoids. Note that the relative average absolute errors between the volumes of the deformable superellipsoids and the segmented data along the cycle are less than 1%.

Figure 12 .

 12 Figure 12. Cylindrical coordinates of two points M and M ′ , which represent the same parametric point at the begining and at the end of the cardiac cycle.

Figure 13 .

 13 Figure 13. Trajectories on a synthetic example. Top: original trajectories. Bottom: estimated trajectories. Left column: view from the apex. Middle column: front view. Right column: view from the base. Dynamic presentation in the video.

Figure 14 .

 14 Figure 14. Twist component estimation on a synthetic example. Left: actual values. Right: estimated values (see text).

Figure 16 .

 16 Figure 16. Amplitude of motion computed at each model point during a cardiac contraction (SPECT). Dynamic presentation in the video.

Figure 15 .

 15 Figure 15. Trajectories of points on the model of the epicardium during the systolic stage (contraction, during three successive time frames) from two viewpoints (SPECT). The two surfaces represent the deformable superellipsoids respectively at the beginning and end of systole. Dynamic presentation in the video.

  for y and z). Equation 11 holds under very general assumptions: u must be continuous in D and have continuous bounded first partial derivatives; the boundary surface ∂ D must have continuously varying tangent planes, except at finitely many vertices and edges. From this equation, we obtain (among others) the following formula: n) d S ,

Table 1 .

 1 Typical computation times with an increasing number of control points of the FFD for 20 iterations.

	Size of the FFD	Computation time (s)

Table 2 .

 2 Averages of errors between original data and parametric models for SPECT data (in mm).

		Separate computation Simultaneous computation Precision loss %
	Epicardium	0.097590	0.102623	5.2
	Endocardium	0.128125	0.135583	5.8

Table 3 .

 3 Tracking the complete DSR and SPECT sequences: number of points (FFD and model), number of iterations, number of time frames and computation times.

		FFD Model Iterations Frames Computation time (min)
	SPECT 125	2500	10	8	25
	DSR	125	5000	10	18	72

Figure 8. A cross-section over time, during a cardiac cycle (DSR). Top: original image. Middle: segmentation. Bottom: reconstruction with a deformable superellipsoid involving 11 parameters + 125 points.

ACKNOWLEDGEMENTS

We would like to thank Grégoire Malandain who contributed to the segmentation of the images presented in this article; and Jérôme Declerck, Serge Benayoun and Alexis Gourdon for constructive remarks; many thanks also to Lewis Griffin and Allen Sanderson for their careful reviews of this article, and to the reviewers of this article for their numerous and relevant remarks. The images of the SPECT sequence were supplied by Prof. Michael Goris, from Stanford University. Thanks to Dr Rich Robb and Dennis P. Hanson of Biomedical Imaging Resource, Mayo Foundation/Clinic for the DSR data. This work was partially supported by Digital Equipment Corporation.


Analyzing 3D Images of the Brain

NICHOLAS AYACHE

OVERVIEW

During the past 5 years, there has been a considerable effort of research in automating the analysis and fusion of multidimensional images, yielding important theoretical and practical new results. These results could now be useful to brain research.

Along these lines, and focusing on 3D images of the brain obtained with CT, MRI, SPECT, and PET techniques, I will present a number of methods useful to produce quantitative measurements necessary for an objective analysis of 3D images of the brain.

Such methods include segmentation, shape analysis, rigid and elastic registration, fusion of multimodal images, analysis of temporal sequences (4D data), modeling, and matching of digital atlases. A large number of references on these topics can be found in Ayache (1995a,b), and specific examples from our research group in [START_REF] Malandain | Improving registration of 3-D medical images using a mechanical based method[END_REF], [START_REF] Subsol | A general scheme for automatically building 3D morphometric anatomical atlases[END_REF][START_REF] Thirion | Fast intensity based non-rigid matching[END_REF].

The following text describes the research tracks to be followed in order to optimize the future exploitation of 3D images of the brain in neuroscience.

NEW VOLUMETRIC IMAGES

It is possible to acquire a variety of 3D images of the brain. They come from different modalities like magnetic resonance imaging, computed tomography imagery, or nuclear medicine imagery. New 3D imagery devices are emerging, like angiographic MRI, magnetoencephalography equipment, or functional MRI.

Volumetric images share the particularity of describing the physical or chemical properties at each point of a studied volume. This information is stored in a discrete 3D matrix I(i, j, k) of voxels (volume elements), called a 3D image. A high resolution MR image of the brain may contain more than 16 million voxels, and this makes its manipulation, visualization, interpretation, and/or comparison with other images or atlases quite difficult by hand.

NEW POTENTIALITIES

In fact, there is a new list of potential tasks which could be performed automatically by some dedicated

NEW MEDICAL IMAGE ANALYSIS, GRAPHICS, AND ROBOTICS ISSUES

To complete these new tasks automatically, it is necessary to solve a number of advanced 3D image analysis research problems (and possibly a number of graphics and/or robotics ones):

3D segmentation of images: the goal is to partition the raw 3D image into regions corresponding to meaningful anatomic structures. It is a prerequisite to most of the tasks listed above. Efficient segmentation requires the modeling and extraction of 3D static or dynamic edges and of 3D texture, as well as the generalization in 3D of digital topology and mathematical morphology.

3D shape representation/characterization: this is mainly a prerequisite to solve the registration and identification needs, but also for efficient visualization. It is necessary to describe nonpolyhedral 3D shapes with a reduced number of intrinsic features. This involves mainly computational and differential geometry.

3D rigid and nonrigid registration and fusion of 3D shapes: once segmented and modeled, new algorithms must be designed to reliably and accurately match such representations together, both in rigid and in nonrigid cases. This is necessary to solve the registration and identification needs. This requires the development and comparison of a number of matching algorithms (tree-search, accumulation, correlation, maximum entropy, etc.), as well as the choice of a class of nonrigid transformations (affine, quadratic, splines, etc.).

3D standard geometric models must be constructed, based on the automatic registration of a number of intersubject anatomic structures or activations. This requires the previous nonrigid matching tools, but also advanced statistical models of random shapes and spatial activations. Once such models are built, it is necessary to build specific registration tools, to adapt these models to the images of an arbitrary subject, allowing the automatic interpretation of the anatomic structures, but also the detection (and quantification) of important deviations from the standards.

Virtual reality, physical models, and image-guided robotics should be developed to provide realistic interactive visualization and to help in planning and simulation, in particular to help in positioning measuring devices within the brain with an extreme accuracy.

CONCLUSION

To exploit fully the potentialities offered by the new 3D images of the brain, there is a strong need for research mainly in the field of 3D image processing and analysis, but also in the fields of 3D graphics and medical robotics.

I believe it is very important that we help, in Europe and in the United States, the best research in these fields, and encourage a strong collaboration between these experts and the experts in neurosciences, to ensure significant progress in the understanding of the brain in the coming decade.

In this spirit, the morphometric analysis of the brain (i.e., the production of quantitative shape measurements on the brain) is the central topic of a European research consortium called Biomorph, including our research group, starting in 1996, and to be funded by the European Community within the Biomed-2 framework.

  software, and which could help brain research significantly: c producing a realistic 3D visualization of the major brain structures in anatomical images; c producing quantitative measurements of shapes and textures in anatomical or functional images; c studying temporal evolution of anatomical or functional images; c fusing of anatomical and functional images; c providing intersubject comparisons of anatomical or functional images; c building geometric models of standard anatomy from anatomical images; c building geometric models of standard activation from functional images; c evaluating deviations from standard anatomy and/or activation; and c planify and accurately position in space a number of physical measuring devices like microelectrodes, for instance.


Andr P Gu

Xavier Pennec

Nicholas Ayache

Medical Image Registration using Geometric Hashing

a pre-processing step to identify various anatomical structures in both images. Such structures are then registered, and the transformation that was found to best superimpose the structures is applied to the underlying images.

We present recent work at INRIA on registering medical images using Geometric Hashing. Here, registering means determining the rigid displacement that aligns two images to obtain the best overlap of anatomical structures represented inside the images.

We analyze two methods that perform hashing on geometric invariants and use a voting procedure, inside a six dimensional accumulator implemented using a one dimensional hash table, to determine the best candidate rigid displacements. The rst method uses mostly unary geometric invariants, with an option to use binary invariants, while the second method uses binary invariants exclusively. The rst method uses characteristic ridge curves extracted on the surface of objects for registration purposes, while the second method uses characteristic extremal points.

We present registration results using both methods, and report experimental registration accuracy data. We also analyze the computational complexity of the registration and compare with the complexity of alternate registration methods. Using data from a joint study with an industrial partner we show that Geometric Hashing can provide a practical, fast, and accurate method for registering CAT scan images of the same patient.

Introduction 1.Registration of Three Dimensional Medical Images

This article develops methods for registering 3D medical images. By registering, we mean determining the rigid displacement (combined rotation and transla-tion) that aligns two images to obtain the best overlap of anatomical structures from the two images. Image registration is recognized as a fundamental problem in Medical Imaging. Fusing and comparing medical images is very useful, and sometimes necessary, for diagnosis, patient follow up, treatment planning and surgical planning. Today a great variety of medical image modalities are available, such as conventional X-ray, X-ray Fluoroscopy, X-ray Computed Tomography (CT), Magnetic Resonance Imagery (MRI), Positron Emission Tomography (PET), Ultrasound (resulting images are sometimes called Sonograms), Single Photon Emission Computed Tomography (SPECT). Without an accurate geometric superimposition, it is hazardous to compare or fuse images corresponding to the same anatomy.

The registration can be performed for a single patient, or across patients (for instance for building a digital anatomical atlas), using a single image modality or multiple modalities. The registration can be rigid, as studied in the present article, or non-rigid, whereby various deformations can be applied to images in order to superimpose them. In general terms, the techniques that have been proposed fall into the following two broad categories: (1) intensity or pixel (voxel) based techniques, and (2) feature based techniques. Intensity based techniques de ne a pixel similarity measure between di erent images ("mutual information", entropy, or simply cross-correlation) and attempt to optimize registration parameters in order to maximize this measure. There is a priori no need to segment the image data, i.e. extract regions, surfaces, curves or points of particular interest. Although the rst intensity based algorithms were generally slow, ine ective at optimizing the registration parameters, and used an objective function of limited relevance, there has been much progress recently in these three respects. Feature based techniques rely on

Our approach

The approach we developed in the Epidaure project at INRIA (France) consists of extracting characteristic curves, called "Crest Lines" on the surface of anatomical structures inside each image, and of registering such lines using Geometric Hashing. In section 2 we develop the method of Gu ziec and Ayache for registering surfaces (and three dimensional images) using characteristic curves 1]. In section 3, we expose the method of Pennec and Ayache 2] that was developed later and builds upon the rst method.

One of Gu ziec and Ayache's original contributions was to exploit curve Fr net frames to directly compute a rotation and translation for each matching pair of curve sample points (instead of using three pairs of matching points as was traditionally done). Also, geometric invariants are computed on curve sample points, and are used for indexing the samples. Such invariants are unary they only need one sample point and its neighborhood to be computed. Instead, Pennec must specify a basis feature with respect to which binary invariants are computed. This increases the computational complexity of the method but makes it more robust. The methods also di er in the way they quantize the invariant space (the accumulator) and determine the contribution of a particular transformation to various accumulator cells. Both presentations are accompanied with experimental results and analyzes of computational complexity.

The main reason for preferring Geometric Hashing over other matching methods is its advantageous computational complexity. This is an important issue in Medical Imaging because of the size and intricacy of structures of interest embedded inside three dimensional images. In section 2.3 we study in detail the complexities of various matching methods for the particular problem of matching curves and motivate our choices: Geometric Hashing is more e cient in general and particularly e cient when one object covers a small portion of the object(s) in the database. Geometric Hashing is very well adapted when one wishes to register a shape against a database of shapes, as the complexity of the algorithm is largely una ected by the size of the database. This capability is used in section 2, as individual curves can be recognized and labeled, as well as in section 3. We believe that the problem of registering against a database is likely to become important in future medical imaging applications using digital anatomical atlases or digital databases of diagnosis images, that are currently being developed as pilot projects.

Registering 3D Medical Images Using Geometric Hashing On Characteristic Curves

This section presents a method for registering 3D medical images using characteristic curves. Details of this method and various experimental results can be consulted in Ayache and Gu ziec's publications 1, 3].

Figure 1: Example of crest lines on the brain surface (front view). Lines in red roughly correspond to sulci whereas blue lines could be interpreted as gyri (Courtesy of G. Subsol).

Extraction of Crest Lines

Crest Lines (also referred to as Ridges or Ridge Curves) are, according to our practical de nition, locations of the surface where the maximum surface principal curvature in absolute value reaches a local maximum in the principal direction of maximum curvature. De nitions can vary slightly, but intuitively, such lines correspond to the surface's most salient features. Sample Crest Lines are shown in Fig. 1. Prior to the late 1980s and early 1990s, the interest in Ridges (Crest Lines) was mostly theoretical, in areas of mathematics related to Catastrophe Theory. Practical applications were then discovered by researchers in Computer Vision, Graphics and Medical Imaging together with the speci cation of algorithms for extracting ridges. Several algorithms were used to extract Crest Lines from medical images. Fully automatic methods were developed by Monga et al and by Thirion and Gourdon. Cutting et al developed a method using a template of Ridges to assist in their extraction from image data. This method was extended by Dean et al (see side-bar).

The approach investigated by Gu ziec was to build a di erentiable representation of the surface using a B-spline surface: by decoupling surface u and v variables, by varying the surface resolution, by imposing constraints on the image gradient, and by using surface regularization parameters compatible with decoupling variables, Gu ziec was able to build surfaces with su cient detail and to generate meaningful Ridges.

Aside from the location of the Ridges, we can com-Di erential Geometry and de nition of Crest Lines pute a series of intrinsic surface properties: principal curvatures and principal directions for a point that is not umbilic. Once Ridges have been extracted and their properties recorded, we can design a very ecient registration algorithm based on the idea of Geometric Hashing. Generally, there has been much more research on the registration of surfaces than on the registration of curves. Although the same techniques could potentially be applied, developing a special purpose algorithm for curves can result in increased eciency and speed. This is what our method does.

Description of the Curve Registration Algorithm

We are given a set of model curves fM i g and a scene curve S. We wish to identify a curve M i which has the largest subset of points in common with S after a rigid transformation (matching); and specify that rigid transformation that minimizes a measure of the deviation between the points that were associated in the two curves (registration). Let m denote the number of model curves, and n the maximum number of sample points in a model or scene curve. Figs. 2 and 3 illustrate our curve matching algorithm. In the preprocessing stage, we process each model curve in turn. For each sample point, we com- Figure 2: Preprocessing: after the computation of three intrinsic parameters (we have used ve for the experiments of Section 2.5), each curve point falls in a cell of an imaginary three dimensional accumulator. The accumulator is implemented with a one dimensional hash table: the mapping from three dimensions to one is performed using a hash function. At recognition time, the same point, although in a di erent position and orientation, will fall in the same accumulator cell.

k 1 k 2 k 1 k 2 k 1 k 2 Model
pute invariant parameters (discussion in section 2.4) and store the sample point inside an accumulator cell according to the parameter values. The principle of the method is the same independently of the number of parameters, or dimensionality of the accumulator. In Fig. 2 three parameters are represented; this should not be construed as meaning that the number of parameters must always be equal to three. In most of our experiments, we have used ve parameters. At this point, it is important to note that accumulators are always implemented using a onedimensional hash table with our approach. In a rst pass, the parameter values are transformed to integers corresponding to a bin number for each dimension (quantization of the invariant space). Rather than constructing a high dimension table that would have space allocated for each potential cell, we use a hash function that maps a tuple of integers (bin numbers) to a single integer code. Typically the function uses a combination of binary shift and xor operations, to attempt covering the full range of the computer representation of an integer (role of the shift) and to distribute entries as much as possible (role of the xor). Such a hash function is represented with a pictograph using language C notation for the binary operators in Fig ( 2). An ideal hash function would distribute entries evenly over the entire integer range.

In the recognition stage, we process each sample point of each scene curve: we compute the set of chosen invariant parameters for each sample point and store (a pointer to) the point inside an accumulator cell corresponding to the parameter values.

Next, pairs of points (M; S) from a scene and a model curves inside the same accumulator cell are tested as potentially corresponding points. A rotation and translation is then computed between the two points and associated frames and entered in a second accumulator, of dimension six. Speci cally, if we are given the Fr net frames (t; n;b) and (t 0 ; n 0 ; b 0 ) corresponding to the points M and S, the rotation R which brings the two frames into correspondence is given simply by the outer product: R = (t 0 ; n 0 ; b 0 )(t; n;b) t :

The translational component is then easily obtained. Although in the particular implementation that is discussed, the quantization of the six dimensional transformation accumulator is determined empirically, we can use statistics or methods for propagating uncertainty measurements to determine an optimal quantization. It is important to note that there is a notion of error zone in the vicinity of a transformation. The error zone is de ned as an ellipsoid, whose six axes lengths correspond to uncertainty measurements. To avoid omitting votes because a particular transformation falls near the boundary of an accumulator cell, if the error zone of that transformation intersects other cells, then a vote is added to

6 parameter transformation: determine hypothetical accumulator cell M i S j stored inside transformation Hash-Table invariant parameters Hash-Table (R,T)
Figure 3: Recognition: each time a point of the scene curves is matched with a model point in the accumulator, the Fr net frames allow to compute the six parameters of a rigid displacement. This displacement falls in an imaginary six dimensional accumulator, which is represented using a second hash table. Buckets of the hash table corresponding to the most populated cells provide hypotheses for the registration parameters. such cells as well. In this implementation, the error zone is the same for each transformation. In the implementation of section 3, the size and shape of uncertainty zones varies with the transformation. For practical purposes, once uncertainty values are determined, the cell width is chosen to be a small multiple of the uncertainty, say three times the uncertainty. This guarantees that a vote is added to at most 2 6 cells for each transformation. It is a good idea to enforce a bound on the number of cells a ected by a given transformation: because of the dimensionality (6) of the accumulator, very large hash tables could be required.

Densely populated cells of this second accumulator are detected as hash table buckets that receive a high number of votes. For each such cell, we compute a re ned rotation and translation, by retrieving the pairs of points that hash to the cell and by determining a new transformation that best registers such points. This can be done by solving a standard least squares problem, or by using robust least squares, or better by using a Kalman lter. The re ned transformation that registers the largest number of points for a pre-speci ed maximum mean squared error is retained. For better stability, it is possible to combine the Geometric Hashing method with methods based on nding closest points. This is essentially what is done in the veri cation step of the algorithm described in the second part of this article.

The complexity of the recognition stage is between O(n) and O(mn 2 ), depending on the level of quantization of the parameter space. The length of a bucket's entry parameter space is at the worst O(mn); the method depends on the existence of parameter variations along the model. The complexity of the pre-processing phase is O(mn). However, this complexity number should be amortized because the hash table that was built can be re-used for a number of registration queries.

Comparison of Approaches to Curve Model Matching

The following approaches have been proposed for registering curves: Schwartz and Sharir's method 4] comprises a linear search step for the best o set between portions of curves. They also provide an innovative and powerful method for curve smoothing with bounded error. Ayache and Faugeras 5] hypothesize and then verify matching hypotheses. Their prediction-veri cation (or alignment) method is reviewed next. Zhang 6] identi es the closest pairs of points between two curves, computes a transformation using such pairs, and recomputes closest point pairs once the transformation has been applied. His method is similar to Besl and McKay's Iterative Closest Point method 7]. Our method is related to the method of Kishon et al 8] who compute the curvature and torsion of the curves and use such parameters for hashing. They use a polygonal representation of the curves, and thus vote for a model and a displacement length, representing a di erence between the arc-length locations of the point s l and the candidate matching point m i;j measured relative to some reference point along each curve's representation. Since our representation of the curves includes a di erentiable structure and thus Fr net frames, we may explicitly compute the rigid transformation as part of the matching process.

We now study the complexity of four methods for curve matching. Our study is related to the comparative study reported in Lamdan the Iterative Closest Point method, uses unary invariant parameter values to select pairs of corresponding points, and computes a rigid transformation for each pair of points using Fr net frames.

In the Alignment method, l signi cant points of the curve S are selected. In the prediction stage, we match each point with a model M i and a point along M i having a similar set of invariant parameters. Each such matching constitutes a hypothesis for which we can compute a rigid transformation. This transformation is used to superimpose the model onto the scene and verify the number of correspondences. One can stop as soon as a su ciently good hypothesis is found or explore all hypotheses and keep the best ones. The computational complexity of the prediction stage is between O(ml log n) and O(mln), depending upon the level of discrimination of the invariant parameters. The complexity of the veri cation stage is O(n). Since the number of hypotheses generated is between O(ml) and O(mln), the overall complexity is between O(mln) and O(mln 2 ).

In the Accumulation Method, for each point on the curve S and every point on a model M i having similar invariant parameters, we compute the rigid transformation that registers them using the corresponding Fr net frames. We register a vote in a quantized sixdimensional rigid transformation accumulator. We then locate the buckets in the accumulator receiving a lot of votes. For each model, the complexity is between O(n log n) and O(n 2 ), depending on the degree of discrimination introduced by the invariant parameters. The di erence with Geometric Hashing is that a binary search is used in parameter space as opposed to hashing into buckets.

In In conclusion, the Geometric Hashing method o ers computational advantages, especially in terms of its sub-linear complexity growth in the number of models.

Choice of Invariant Parameters

We rst considered curvature and torsion as only invariant parameters. Curvature and torsion measured on a smooth curve are invariant to rigid displacements. Then the method was extended to semidi erential invariants, combining distance and angles measured with respect to a reference point. Statistical studies indicated that such semi-di erential invariants could be more stable than curvature and torsion alone. This approach was later exploited by Pajdla and van Gool 10].

A few points on the curves called basis points must be selected and used to compute the parameters of the other curve points. For instance, we can compute angles, a Euclidean distance or an arc-length measure between any given curve point and a basis point. Choosing such basis points can be di cult; an interesting possibility is to use maxima of curvature. Also, since the parameters that are used for recognition involve pairs of points, the complexity of the algorithm is increased, and is in principle O(n 2 ) or higher, where n is the number of points in the curves.

Exploiting the property that our curves are embedded on surfaces, we discovered additional intrinsic parameters, di erent than curvature and torsion. The algorithm extracting Ridges computes the surface normal N, the maximum and minimum surface principal curvatures k 1 and k 2 as well as the associated directions e 1 and e 2 (or their opposite). For each point of a curve that we consider, two frames are available (Fig. 4), (1) a surface frame (e 1 ; e 2 ; N) and a curve Fr net frame (t; n;b) where t is the curve tangent, n the curve normal and b the binormal. A review of these notions can be found in a good textbook on Di erential Geometry. We expect that n and N will be di erent as well as t and e 1 . Otherwise, the ridge would be a line of curvature, which is not expected in general. For each curve point we can compute the Fr net frame (t; n;b) where t is the curve tangent, n the curve normal and b the binormal and the local surface frame (e 1 ; e 2 ; N), where e 1 and e 2 are the principal directions of curvature in the surface tangent plane and N is the surface normal.

Our assumption is that we can choose curve sample points that are di erent from surface umbilic points, where e 1 and e 2 would not be de ned, and where k 1 and k 2 would be equal. Various angles and parameters can be computed using both frames, the curve curvature and torsion and the surface principal curvatures. Not all the parameters that can be extracted are independent. We have identi ed a few combinations of ve parameters. One of such combinations of parameters is the following: (1) curvature of the curve k, (2) torsion of the curve , (3) geodesic torsion of the curve with respect to the surface g , and (4,5) two angles measured between the surface frame and the curve Fr net frame:

( = ( d n;N) = ( d t;e 1 )
The two invariants k and g can be substituted with the principal curvatures k 1 and k 2 using the relations: k cos = k 1 cos 2 + k 2 sin 2 g = (k 2 ? k 1 ) cos sin thus providing other combinations of invariants. This raises the interesting question of whether the invariant parameters that we use are indeed independent: we observe that several surfaces can be locally di erent and have locally the exact same ridge (excluding obvious global di erences). A very simple example is given by an ellipsoid E 1 of axes lengths (a > b > c > 0) and an ellipsoid E 2 of axes lengths (a > b > c=2 > 0). The longest (and sharpest) ridge will be precisely the same for E 1 and E 2 but the underlying ellipsoids will have quite di erent principal curvatures. Armed with this observation, we argue that in general the surface frame (e 1 ; e 2 ; N) and the ridge Fr net frame (t; n;b) should not be related, and that the curvature and torsion of the curve should not explain the surface curvature entirely. Consider a surface ridge on an arbitrary smooth surface: while keeping the ridge in its exact same position, by analogy with the ellipsoid example, we can start pinching the surface and changing the local curvature and principal directions for various locations on the ridge. Since the curvature and torsion of the curve alone cannot explain the local surface geometry, angles relating the surface and curve frames must be specied, as well as the portion of the surface curvature not explained by the curve curvature, providing ve independent invariant parameters.

In all cases, the curve Fr net frame must be computed independently from the algorithm that extracts ridges, by di erentiating the curve. In order to compute the di erential parameters of the curve, we build a di erentiable representation using B-splines. This is a standard procedure for which we have added a few improvements. In particular, we have constrained the curve to be embedded in the surface by minimizing the sum of square products between curve tangent and surface normal. An additional bene t of this method is to apply a certain amount of smoothing to the parameters k 1 , k 2 and g as well as to the Ridge position, thereby removing some of the noise. The details of the curve tting method can be consulted in 1].

Experimental Results

We consider two CT scan images of a skull phantom at a resolution of 256 by 256. The skull phantom was scanned in two di erent positions and orientations. Such images were processed by the algorithm of Thirion and Gourdon, resulting in the the sets of ridges displayed in Fig. 5.a. The algorithm described in this section was then used to register the lines (Fig. 5

.b).

We performed registration experiments using two CT images of a vertebra scanned in two di erent positions and orientations. The result of ridge extraction is shown in Fig. 6.a. The result of registering the ridges using our algorithm is shown in Fig. 6.b.

Within the framework of a joint study with G.E. Medical Systems France, our goal was to assess registration accuracy for the particular application of single modality, single patient registration. This partic- ular problem has several important real world applications. This technology applies for post-operative or post-treatment assessment and patient follow-up studies. Another potential application consists of fusing data-sets acquired in di erent orientations, to compensate for the anisotropy of the process of building 3D volume data by stacking scanned 2D slices.

We used a calibration phantom speci cally designed for the study. The phantom was built using a skull phantom and an a xed pair of N shaped tubular markers. The phantom was imaged (CT scanned) in two di erent positions and orientations (256 x 256

x 30 images with 1mm slices). The second position was registered to the rst position using our method of ridge curve matching. We then computed the distance between points sampled on the N shaped markers in the rst view and corresponding points on the second view after registration. Such distances were found to be between 0.05 mm at the N shaped marker centers and 0.85 mm at the marker peripheries. Most of the areas of clinical interest were situated in the vicinity of the N shaped marker centers, indicating that the 0.05 mm registration error would more likely apply. The registration technology that was described in the present article was transferred to the company Aleph Medical in Grenoble (currently: Focus Medical).

Substructure Matching with Geometric Hashing on Frame Features

In the previous section, we have rst modeled objects (anatomical structures) using surfaces and then using curves embedded on the surfaces. In order to simplify the matching and registration process, it is often necessary to go one step further and extract isolated landmarks. Indeed, if we could select on each line a mean number of l features that are supposed to be well conserved, we could focus on these l features only instead of the n points of the whole line. The complexity of the previous algorithm would then be reduced from a lower bound of O(n) to O(l). Moreover, isolated features, or landmarks, could appear even in the absence of curves, as with the modeling of proteins in section 3.2. However, the increase in the number of lines from the previous section and the increase of the noise with real patient images drastically augments the density of invariants in the hash table. This could lead to an important number of false positives that overwhelm the correct matches. The maximum complexity would then be reached and the algorithm could even provide a wrong answer.

In order to have a correct algorithm, we show in this section how to keep a highly selective hash table by using binary invariants instead of just unary ones. This approach is related to the semi-di erential invariants methods previously discussed. The drawback is an increase of the complexity from O(n) to O(l 2 ) (in practice, we observe that the lower bound is respected and the correct solution is always found).

Registering Medical Images with Frames

To obtain precisely located features, Thirion introduced new feature points called extremal points, which are points of locally maximal curvature in both principal directions of the surface. These points are speci c points on the crest lines and can be seen as the generalization of corner points for smooth surfaces in 3D. Another way to view such points is to consider extremal lines instead of just crest lines. Extremal lines are de ned as the loci on the surface where at least one of the principal curvatures is extremum in the associated principal direction. Extremal points are at the intersection of these lines. To the position of the point, we can add the normal vector N and the two principal directions e 1 and e 2 of the surface to constitute a local coordinate system or a frame (Fig. 7). In this context, each medical image is modeled by a set of frames and the matching problem is to nd correspondences between a subset of these frame sets in di erent images. We have developed in 11] methods to solve the dual problem of registration of matched frames (i.e. nding the rigid transformation).

Preprocessing

Step To obtain an invariant representation with respect to the global position and orientation of the considered structure, we can express the con guration of all frames relative to one frame (the basis). For e ciency, we store this representation in a hash table and, for correctness, we include the uncertainty of each invariant. As only part of the frames are in the same con guration in the two images, the one chosen as the basis may not be present in the other image. The preprocessing step is thus repeated with each frame as the basis.

Recognition

Step We choose a frame of the second structure (the scene) to be the basis for computing the invariant representation and retrieve, thanks to the hash table, what are the compatible model frame couples. If the basis belongs to a common substructure, then a signi cant number of frames are in the same con guration with respect to it. We then match the model and scene bases (Fig. 8). This process is repeated for every amino acid as the basis to nd its possible matches in the model and we only keep the matches that are above a given threshold (typically 4 or 5 in our applications).

Clustering Compatible Matches and Veri cation. For each individual match, we maintain dur- Indexing the invariant space: Hash table Figure 8: Preprocessing: the 6D invariant vector associated with every couple of model frames is computed with its error zone and used as an index for the couple in the hash table. Recognition: for each scene frame couple, we compute the 6D invariant vector and retrieve through the hash table every compatible model frame couple. For each such couple, we tally a vote for the matching of the reference frames (here the match (F m i , Fs j ) scores 2).

ing the recognition step an estimation of the associated transformation by fusing the transformations between con rming frames. To group matches belonging to the same (rigid) substructure, we run a very simple clustering algorithm on the associated transformation. Indeed, up to measurement errors, frames should undergo a similar transformation within a single substructure. Each cluster is then re ned by an iterative closest neighbor technique where we enforce symmetry of the matches and verify their validity with a 2 test.

In the case of medical images, matching model and scene crest lines is not su cient since di erent line matches can correspond to di erent transformations. Thus, we run once again our clustering algorithm on transformations to nd out the compatible line matches and we obtain a single transformation from the model to the scene image.

Registration of CT images of the skull We present in Fig. 9 an example of the registration of two CT images of the skull of the same patient acquired in di erent positions. About 75 crest lines are matched with more than 4 extremal points among the 550 in each image leading to a total of 550 matched extremal points (only on the 75 matched lines). Using the techniques described in 12], we have computed that the registration accuracy (the expected standard RMS on image super-imposition in the area of matched features) is 0.04 mm.

Registration of MR images of the head Fig. 10 is an example of registration of two MR images of the same patient. In this case, 80 crest lines are matched among the (about) 620 in each image, for a total of 350 matched extremal points. The accuracy of the registration is 0.06 mm. In this case, we have 24 images of the same patient, which allow us to run the statistical validation method we developed in 12]. It con rms our accuracy estimation.

Application to Molecular Biology

The recognition algorithm presented is this section was originally developed to nd common substructures between proteins 2]. We show very brie y how this problem can be put in the same framework and show an application example.

Non-related proteins can have very di erent structures, but when they share a common biological function, often they share a common 3D substructure that realizes this function: this is a motif. Binding motifs can be widely distributed in the protein sequence and involve only few, if any, consecutive amino acids. Thus we cannot really use the sequence and we have to focus on the pure 3D structure. Moreover, the nature of amino acids themselves can vary since only part of the residues are used to bind. The only meaningful information is thus the geometric con guration of the amino acids in space. This observation leads Fischer et al 13] to apply the geometric hashing technique to this problem.

Each amino acid has 4 atoms participating to the backbone of the protein, 3 of them always being in the same geometric con guration. We use these 3 atoms to de ne the con guration of the amino acid in space (Fig. 12). Each amino acid is thus modeled as a frame (as were extremal points in medical images), and a protein is modeled by an unordered set of frames (as were crest lines before). To nd similar substructures and 2WRP: the HTH motif. We can see that not only the backbone is very well matched, but also collateral chains are pretty well conserved.

in two proteins, we now have to nd two subsets of frames that are in the same con guration, up to a global rigid transformation. Since this is exactly the problem we address previously for medical images, we can use exactly the same algorithm. We next study two proteins known to bind to the DNA in two di erent organisms: the tryptophan repressor of E. Coli and CRO protein of phage 434. These two proteins bind tightly to the DNA at a gene's promoter, thus preventing RNA polymerase from xing and hence blocking the initiation of the transcription (they are repressors). They should thus have a common substructure that realizes this biological function.

The comparison of the two previous proteins gives one principal common substructure which is, as expected, the Helix-Turn-Helix motif, responsible for the binding to the DNA. The registration of the HTH motifs in the two proteins shows that not only the backbone is very well matched, but collateral chains are also well conserved. The use of frames instead of just points allows the algorithm to discard several false matches of amino acids (4 in the HTH case). Moreover, it makes the algorithm robust with respect to parameters that tune the quality of the matches (uncertainty attached to frames on input).

Discussion

In the rst part, we have build upon previous methods for registering curves to incorporate geometric invariants associated to a ribbon, or strip, of surface. Following the indexing step of Geometric Hashing in the space of invariants, we have introduced a six dimensional accumulator to e ciently select the rigid transformation that register the maximum number of curve points. Our algorithm thus acts as a combination of two ltering procedures for isolating the best candidate transformations. The computational complexity is somewhere between O(n) and O(n 2 ), where n is the number of sample points chosen on the curves to be registered. Experimental results using data from a joint study with a Medical Imaging manufacturer indicate that the method is suitable for registering same-patient CAT scans because of its repeatability, speed, and accuracy. In the second part, we have developed a new extension of geometric hashing based on frame features. We show that Geometric Hashing can be successfully applied to very di erent problems such as substructure matching of proteins and matching of medical images.

The examples shown in this article indicate that Geometric Hashing can provide a powerful method for registering highly detailed three dimensional medical images of the same modality (CT registered against CT or MR against MR). The levels of registration accuracy that can be obtained are compatible with clinical requirements; also, the computational performance is very good, which allows to perform numerous registration attempts and favors experimentation. Aside from registration, Geometric Hashing also works for recognition: we believe that Geometric Hashing is a promising method for registering against a database of images, and for automati-cally recognizing a feature, such as a speci c anatomic structure, or the absence thereof, inside an image. Also, Geometric Hashing provides useful statistics on the number of features that were paired and on how well they could be registered; this information, which is unbiased by a user-speci ed initialization, can be used for evaluating the registration.

Our experience with the method has taught us a few limitations: a primary requirement is to isolate corresponding geometric features: low level of detail imaging modalities, such as SPECT and PET lack geometric features that Geometric Hashing uses; generally, it is di cult to extract reproducible features for multi-modality registration; for the case of deformable registration, the di culty of nding invariant parameters hampers the method. As discussed in Section 3, when there are few corresponding features, or when there is signi cant noise corrupting the features, the method can break down unless serious uncertainty handling techniques are implemented.

Figure 4 :

 4 Figure 4: Surface ribbon centered around a space curve. For each curve point we can compute the Fr net frame

Figure 5 :

 5 Figure 5: (a) Left: Crest lines on two 3D CT images of a skull obtained with the algorithm of Thirion and Gourdon. (b) Right: Superimposition of the lines of a. with our software in 7.7 seconds CPU on a DEC 5900.8 curves are matched out of 15 and 32 curves. 53 points are matched out of 666 and 1124. Note that these experiments were performed before the DEC alpha architecture was available. CPU timings would be signi cantly decreased today.

Figure 6 :

 6 Figure 6: (a) Left: Crest lines on two 3D CT images of a vertebra. (b) Right: Superimposition of the lines of a.with our software in 8.9 seconds CPU on a DEC 5900. 12 curves were identi ed out of 28 and 33. 82 points were identi ed among these curves out of 1250 and 1213 points respectively in the two di erent sets of curves.

Figure 7 :

 7 Figure 7: Crest line with an extremal point, modeled using a frame.

Figure 9 :

 9 Figure 9: Example of registered crest lines between two CT skull images of the same patient. Extremal points are represented by a color change from yellow to blue on the lines. (a) Left: Front view with all crest lines from the two skulls. (b) Middle: Left view of the matched crest lines. (c) Right: Closeup on the occipital foramen on the right. In this last image, the width of a line is a tenth of a voxel, which shows the very precise registration of these extremal points. One can also see that the trihedron part of the matched frames is very well conserved.

Figure 10 :

 10 Figure 10: Example of registered crest lines between two MR images of the brain (same patient). Extremal points are represented by a color change from yellow to blue on the lines. (a) Left: View of matched crest lines from the top. (b) Middle: view from the right of the head. (c) Right: closeup to show the quality of the frame matches.

Figure 11 :

 11 Figure 11: (a) Left: The CRO protein (2CRO) of phage 434. (b) Middle: The tryptophan repressor of E. Coli (2WRP). The matched part is displayed in yellow. (c) Right: Registered matches found between 2CRO

Figure 12 :

 12 Figure 12: Modeling an amino acid by a frame.

  .P. do Carmo. Di erential Geometry of Curves and Surface. Prentice Hall, Englewood Cli s, New Jersey, 1976. I. Porteous. Geometric Di erentiation: For the Intelligence of Curves and Surfaces. Cambridge University Press, 1994. J.J. Koenderink. Solid Shape. M.I.T. Press, 1990.

	Extraction of Crest Lines and Extremal Points
	O. Monga et al. Using third order derivatives to extract ridge lines in 3-D images. In Proc. CVPR'92, Urbana Champain (Il), June 1992.
	J-P. Thirion and A. Gourdon. The 3-D marching lines algorithm. Graphical Models and Image Processing, 58(6):503 509, 1996.
	C.B. Cutting. Applications of computer graphics to the evaluation and treatment of major craniofacial malfor-mation. In Udupa and Herman, editors, 3-D Imaging in Medicine. CRC Press, 1989.
	A. Gu ziec. Surface representation with deformable splines: Using decoupled variables. IEEE Computational Science and Engineering Magazine, 2(1):69 80, 1995.
	J-P Thirion. New feature points based on geometric in-variants for 3-D image registration. Int. J. Computer Vision, 18(2):121 137, 1996.

M

Table 1 :

 1 Lower and upper complexity bounds for the four di erent methods of curve recognition: Alignment, Accumulation, Geometric Hashing and ICP.

	and Wolfson 9], with the following di erences: each method, except

  the Iterative Closest Point method (ICP), for each n points on the scene curve, the closest point on a particular model is determined in O(log n) time in average and in O(n) in the worst case. Once closest

	points are determined, a transformation is computed in O(n) that registers the corresponding points. The entire operation is then repeated a suitable number of times i. The complexity is between O(mn log n) and O(imn 2 ). It is also proportional to the number of initial positions.

The complexity results are summarized in table 1.

Is Scene Base Is j,1 Is j,2 Fm 2 Im i,3 Fm 1 Fm 4 Fm i Fm 3 Model Base Im i,1 Im i,4 Im i,2

			Fs 1
		j,3	Fs j
	Fs	3
			Fs 2

Acknowledgements

We wish to thank J. Knoplioch of General Electric Medical Systems Europe, for providing data sets used in studying the accuracy of registration, as well as J.L. Coatrieux, for providing CT images of a vertebra.