Serge Benayoun

Nicholas Ayache

Dense Non-Rigid Motion Estimation in Sequences of Medical Images Using Differential Constraints

Keywords: medical imaging, adaptive meshes, non-rigid motion, differential geometry, finite element method

We describe a new method for computing the displacement vector field in time sequences of 2D or 3D images (4D data). The method is energy-minimizing on the space of correspondence functions; the energy is split into two terms, with one term matching differential singularities in the images, and the other constraining the regularity of the field. In order to reduce the computational time of the motion estimation, we use an adaptive image mesh, the resolution of which depends on the value of the gradient intensity. We solve numerically the minimization problem with the finite element method which gives a continuous approximation of the solution. We present experimental results on synthetic data and on medical images and we show how to use these results for analyzing cardiac deformations.

Introduction

The problem of computing non-rigid motion from digital images is one of the most attractive challenges in computer vision. There are many applications including, for example, model-based image compression [START_REF] Tang | Quantifying facial expressions: smiles[END_REF] or medical diagnosis [START_REF] Ayache | Medical computer vision, virtual reality and robotics[END_REF]. Medical imaging produces sequences of 2D and 3D images which describe the time evolution of anatomic structures. This paper particularly deals with heart images. Cardiovascular diseases represent the "greatest scourge afflicting the population of the industrialized nations" (Braunwald, 1992). Medical imaging now produces a lot of modalities in order to visualize the cardiac cycle, like Magnetic Resonance Imagery (MRI), Nuclear Medicine Imagery (NMI) or Computed Tomography (CT). Because the heart is a nonrigid organ, it is difficult to assess its deformations. This objective, however, seems very important for a better understanding of this organ. For example, some pathologies characterized by ischemic or infarcted tissues are detectable because of an abnormal local motion. Image Processing enables an accurate computation of the motion of the cavities, which could help cardiologists to localize such tissues and/or to evaluate the effects of a medicinal treatment.

Many works have already been done on this subject. The most commonly used method consists of recovering the motion by tracking a deformable model [START_REF] Terzopoulos | Constraints on deformable models: recovering 3D shape and nonrigid motion[END_REF][START_REF] Cohen | Tracking points on deformable objects using curvature information[END_REF][START_REF] Nastar | Vibration modes for nonrigid motion analysis in 3D images[END_REF][START_REF] Park | Volumetric deformable models with parameter functions: a new approach to the 3D motion analysis of the LV from MRI-SPAMM[END_REF][START_REF] Bardinet | Superquadrics and free-form deformations: a global model to fit and track 3D medical data[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]. It needs a twostage approach. First, one fits contour points with a deformable model. It can be either a free shape model like snakes [START_REF] Kass | Snakes: Active contour models[END_REF][START_REF] Cohen | Finite element methods for active contour models and balloons for 2-D and 3-D images[END_REF] or a constrained model like superquadrics [START_REF] Terzopoulos | Dynamic 3D models with local and global deformations: Deformable superquadrics[END_REF][START_REF] Bardinet | Superquadrics and free-form deformations: a global model to fit and track 3D medical data[END_REF]. In this latter case it is necessary to add local deformations in order to recover complex shapes. However, free shape models uses a smoothness constraint for controling the regularity of the fit. The second stage consists in tracking the model over the time sequence. It is based on the matching of some particular points based on metric or geometric similarity measures. Our method mainly differs from these works by the fact we do not use any shape model. Several reasons explain this choice. First it can be very difficult to recover anatomic organ which may have complex shape and topology. Second, the motion information may be not only present in one structure but everywhere. The last reason concerns deformable models which present some specific problems, one of them being that they have to be initialized very close to the solution.

In this paper, we present a method which computes the displacement vector field between two successive frames without prior segmentation. This procedure is applied sequentially to the set of frames. We pursue the original idea introduced in a previous paper [START_REF] Cohen | Tracking points on deformable objects using curvature information[END_REF] of using geometric singularities to anchor the motion field on a few but reliable characteristic points, and to propagate this sparse motion field, first to contours, and then to every image point. The idea of using salient features to guide the components of the nonrigid motion field was also suggested in different studies reported by [START_REF] Amini | Bending and stretching models for LV wall motion analysis from curves and surfaces[END_REF] and [START_REF] Benayoun | Local quantitative measurements for cardiac motion analysis[END_REF]. The implementation of our method is based on an energy minimization. We first consider a mesh whose resolution depends on the value of the gradient intensity computed on the data [START_REF] Terzopoulos | Sampling and reconstruction with adaptive meshes[END_REF]. This adaptive mesh mainly allows to reduce the computational requirements without decreasing the accuracy of the displacement field.

The paper is organized as follows. First we describe the general formulation of our problem (Section 2). Next we present our model of adaptive image meshing (Section 3). The tracking stage is next modeled as an energy minimization problem (Section 4). It is solved by a mathematically rigorous finite element method. Finally we show experimental results on synthetic data, and on time sequences of 2D and 3D medical images (Section 6). We make use of analytical modal analysis [START_REF] Nastar | Vibration modes for nonrigid motion analysis in 3D images[END_REF] as an instantaneous method for analysis and smoothing of the computed displacement field, and also show such results in this experimental section.

Problem Formulation

Our aim is to compute the displacement vector field between two successive frames of a sequence without prior segmentation. This means that we want to determine for every material point in a frame its new position in the following frame. The fundamental hypothesis we apply is that the mean displacement amplitude between successive frames is low. Under this constraint, it becomes possible to track a point just by searching for a similar point in the following frame. Two points are said similar if they share some features computed from the gray level data. The most straightforward choice consists of using the image intensity. However, in general, the hypothesis of a perfect conservation of the intensity is not robust. A most robust approach assumes the conservation of contours. Contour points are generally detected by extracting the local maxima of the intensity gradient amplitude in the gradient direction [START_REF] Monga | Recursive filtering and edge closing: two primary tools for 3D edge detection[END_REF]. They are also characterized by a high intensity gradient amplitude. But matching points is not sufficient neither to recover the tangential component of the motion along the contours, nor to compute a dense motion field.

In order to recover the tangential component of the motion, we look for high curvature points along contours. For example, the points belonging to ridge lines on surfaces correspond to the local maxima of the principal curvature in the associated principal curvature direction. They are anatomically significant in medical imaging and stable under small deformations [START_REF] Monga | Using partial derivatives of 3D images to extract typical surface features[END_REF]. Thus they represent good landmarks for tracking.

A dense motion field is finally recovered by imposing a regularity constraint everywhere, which provides a propagation of the tangential component of the motion from high curvature points along the rest of the contours, and a propagation of the motion field from contours to the overall image.

To reduce the computational requirements, and because we do not need a very dense solution far from the contours, we introduce an adaptive image mesh. This mesh has a structure similar to the original frame, and keeps its topological and geometrical properties. But the number of nodes is reduced relatively to the original number of pixels or voxels. The density of nodes is increased near the contours where the underlying computed field is meaningful. The problem now reduces itself to searching for every node of the mesh its new position in the following image.

Having defined our constraints, we can now formulate the matching process as a minimization problem. Let f be the correspondence function which matches nodes in one frame with points of the following frame; then we define an energy E e ( f ) which measures the quality of the match. Because the problem of minimizing this energy may have locally several equivalent solutions, we add a second energy E i ( f ) which constrains the regularity of the function f . It means that, when it is possible to match a node with several points, the most regular solution is favored. We call E e ( f ) the external energy because it depends on image information, and E i ( f ) the internal energy because it concerns only internal properties of f .

Adaptive Meshes

The use of adaptive meshes is a classical tool in computational fluid dynamics. They have already been used in computer vision, specifically for range data reconstruction [START_REF] Terzopoulos | Sampling and reconstruction with adaptive meshes[END_REF] or surface mesh decimation [START_REF] Gourdon | Simplification of irregular surfaces meshes in 3D medical images[END_REF]. One of the advantages of this technique is to decrease the numerical complexity of a structure while preserving its accuracy. In our case, it will be mainly used to reduce the original frame structure. In particular, for 3D images, it would be too long to compute the displacement of every voxel.

A mesh is a set of nodes connected by springs. Each spring models the connection between two neighboring nodes by the relationship:

S i j = k i j d i j X i j (1)
The strength vector S i j exerted by the node j on node i is in the direction of the unit vector X i j which connects the two nodes. The intensity is proportional to the elongation of the spring d i j , and also to its stiffness k i j .

There are different techniques for building up an adaptive mesh. One generally makes the distinction between unstructured mesh and structured mesh. For unstructured meshes, the adjacency graph of the nodes is not invariant. For example, a 3D mesh is initialized at the higher resolution (voxel resolution) and irrelevant nodes are progressively eliminated. This technique is suited for example to data compression because no further computation follows the mesh adaptation. Since it is clearly not our case, we prefer to use structured meshes where connectivity between nodes is preserved. For 2D images, we use a rectangular mesh with a 8-neighborhood system. For 3D images, we use a 18-neighborhood system. At the beginning of the process, the mesh is regular and next the nodes move towards relevant points of the image (see Figs. 1 and2).

Our aim is to increase the density of nodes near the high intensity gradient points of the image. This is because of our motion computation method which is well defined at these feature points. One way of doing this is to introduce an image dependent force. This force pulls the nodes towards high intensity gradient points.

It is derived from a potential field computed over the whole image. In [START_REF] Terzopoulos | Sampling and reconstruction with adaptive meshes[END_REF], the adaptation is due to the spring stiffness which is made dependent on the gray level gradient. We also implemented this technique and results are very similar.

More precisely the adaptation of the mesh results from the solution of the system of dynamic equations. Each equation models the spatial behavior of a node when it is influenced by an external force:

m i d 2 P i dt 2 + γ i d P i dt = F i , i = 1 . . . n (2)
where n is the number of nodes, P i ∈ IR 3 is the location of the node i, m i its mass, γ i its damping coefficient, and F i ∈ IR 3 the sum of the applied forces. Notice that we impose an infinite mass for nodes which are in the border of the mesh. Thus they cannot move and we are sure that the mesh covers completely the image. F i is composed of two terms. First, a spring term S i = j∈V i S i j where V i is the 18-neighborhood of node i (see Fig. 3). Second, an image term G i which attracts the node towards high intensity gradient points.

Computation of G i .

Using an edge detector we obtain a binary image with 1 at contour points locations and 0 elsewhere [START_REF] Monga | Recursive filtering and edge closing: two primary tools for 3D edge detection[END_REF]. Next we apply a thinning operation to the edges so that they are all a single pixel wide. Finally, we smooth this binary image and we compute G i by differentiation of the smoothed image. Because we use a convolution filter for differentiating (see Appendix A.1), the direction of the force is a combination of many directions but with a predominance of the direction(s) towards closest contour point(s). A similar result would have been obtained by directly differentiating the grey level function. However, it has the disadvantage that the amplitude of the force may vary along a contour because of its dependence on the intensity gradient.

We use a finite difference scheme to solve the system of Eqs. (2). The algorithm is the same explicit Euler method used in [START_REF] Terzopoulos | Sampling and reconstruction with adaptive meshes[END_REF]. In our experiments, we choose to consider a node each four or five voxels. The mesh deforms through application of the image forces and converges to an equilibrium state.

Energy Minimization

Consider two successive images in a time sequence. We denote by 1 (resp. 2 ) the domain of the first (resp. second) frame. The correspondence function f is a vector function which associates each 3D point

(x 1 , x 2 , x 3 ) ∈ 1 with another 3D point (X 1 , X 2 , X 3 ) = f (x 1 , x 2 , x 3 ) ∈ 2 .

External Energy

The correspondences are achieved by minimizing an external energy E e which preserves differential singularities computed on the images. This energy is split into two terms, a curvature dependent energy E c and a gradient dependent energy E g :

E c ( f ) = 1 (C 1 (x) -C 2 ( f (x))) 2 R c (x) dx = 1 e 2 c (x)R c (x) dx
where x = (x 1 , x 2 , x 3 ) and C 1 (resp. C 2 ) is a curvature dependent function computed in the first (resp. second) image and R c is a curvature dependent weighting function. Therefore, e c is an energy which will measure the discrepancy between contour curvatures at points where R c is different from zero.

E g ( f ) = 1 (N 1 (x) -N 2 ( f (x))) 2 R g (x) dx = 1 e 2 g (x)R g (x) dx
where N 1 (resp. N 2 ) is the gradient amplitude dependent function computed in the first (resp. second) image, and R g is a gradient dependent weighting function. Therefore, e g is an energy which will measure the discrepancy between the gradient amplitudes along contour points when R g is different from zero. The mathematical formula of the curvature dependent functions are given in the Appendixes A.3 and A.4. Concerning the gradient dependent functions, the formula can be found in [START_REF] Monga | Recursive filtering and edge closing: two primary tools for 3D edge detection[END_REF].

Internal Energy

The problem of minimizing this energy is an ill-posed problem because there are many possible solutions. So, we add an internal energy E i which constrains the regularity of the function:

E i ( f ) = 1 ∇ f (x) 2 dx (3)
This also allows the displacement field to be propagated in image areas where the previous energies are not relevant. Finally we solve the minimization problem:

Min f E c ( f ) + E g ( f ) + E i ( f ) (4)

Weighting Functions

The internal and external energies are balanced through the two weighting functions R g and R c . R g (resp. R c ) is a normalized increasing function of the gradient amplitude (resp. curvature):

R g (x) = h(N 1 (x)) R c (x) = h(C 1 (x))
where h acts like a high pass filter. Up to now we used the following form for h:

h(v) = -2 v -v m v M -v m 3 + 3 v -v m v M -v m 2 , (5)
where v m and v M are the minimal and maximal values. This choice is due to the following considerations.

Weighting functions must locally enforce the external gradient and curvature energies at points where gradient and curvature values are significant. Conversely, they inhibit these external energies when a point has no significant gradient or curvature features (see Fig. 4).

In that case, only the internal energy constrains the solution.

In fact, we have chosen to convert R g into R g :

R g (x) = (1 -R c (x))R g (x)
in order to make the curvature energy dominate the gradient energy. This combination is justified for high curvature points because it favors the attraction by another high curvature point. Its effect is illustrated in Fig. 5.

Partial Differential Equation

A necessary condition for the optimality of f = (u, v, w) is the Euler-Lagrange equation ∇ E( f ) = 0 which is equivalent to the partial differential equations:

     -u = R g ∂ N 2 ∂ x (u, v, w)e g + R c ∂C 2 ∂x (u, v, w)e c -v = R g ∂ N 2 ∂ y (u, v, w)e g + R c ∂C 2 ∂y (u, v, w)e c -w = R g ∂ N 2 ∂z (u, v, w)e g + R c ∂C 2 ∂z (u, v, w)e c with boundary conditions f | ∂ 1 = f b .
Consider the first equation; as its second member also depends on the solution, we study the associated evolution problem:

∂u ∂t -u = R g ∂ N 2 ∂ x (D)e g + R c ∂C 2 ∂ x (D)e c (6) with    D = (u, v, w) (u 0 , v 0 , w 0 ) = (u, v, w) t=0 f | ∂ 1 = f b
A solution of the original equation will correspond to a stationary solution of the evolution problem. Moreover, the introduction of a temporal variable allows us to reconsider the using of weighting functions. It becomes possible in particular to reduce little by little the influence of the curvature constraint. One way to do this is to make R c dependent on time:

R c (x, t) = R c (x) 1 1 + αt
where α controls the variation speed.

Indeed for high local deformations, the curvature is not invariant although it remains useful to guide the matching of singular points. At the end of the convergence process (6), the elimination of the curvature constraint ensures that all contour points are well matched (see Fig. 6). The numerical solution of the Eq. ( 6) can be found in [START_REF] Ciarlet | The Finite Element Methods for Elliptic Probems[END_REF].

Experimental Results

We present experiments on sequences of 2D and 3D medical image.

A Synthetic Example

We first illustrate our algorithm by results on synthetic data. We build a 3D image representing a cube. The size of the image is 64 × 64 × 64. By application of an affine transformation, we obtain a second image (the mean displacement issued from this transformation is almost 3.5 voxels). We next consider a 15 × 15 × 15 mesh. After the adaptation stage, the nodes are concentrated near the high gradient points of the cube image. Next, we compute correspondences between the nodes and the points in the deformed cube image. In order to check the accuracy of the tracking, we then compute for each node the distance between the displacement obtained with our algorithm and the exact displacement given by the affine transformation. The results are presented on Table 1. The mean distance is about one half voxel. Notice that the accuracy is better for high curvature nodes, for which the standard deviation of the distance is the lowest.

An Articulated Structure: The Knee

We applied our algorithm to a knee sequence. Sagittal cross-sections have been obtained with a MRI device. 1 As the quality of these images is good, contour points may be well detected. However, the mean amplitude of the displacement between two successive images is quite high for a local approach (about ten pixels). The size of the images is 256 × 256 and we consider a mesh with 46 × 46 nodes. It means that the computational time of the motion field estimation stage is approximately reduced by a factor 25.

Figure 7 shows two successive images on this sequence, the mesh adapted to the first image and the motion field computed between them.

We used mapping to estimate the accuracy of our motion fields. We will notice i ⇐ j the mapping after resampling of the jth image in a sequence onto the ith image. This mapping is based on solving for the mapping f using our proposed method. The mapping algorithm is as follows:

• Consider each pixel p on the ith image • By application of the successive motion fields obtain a pixel q on the jth image (since the motion is only computed on nodes, the displacement on another location is obtained by linear interpolation) • Set the gray level on q to p By this way we obtain a new image which has to be very similar to the original one. The accuracy of the mapping may be checked by computing the difference between these two images. However, it would suppose that the gray level is invariant, which is not true. So we prefer to overlay into the resulting image the contours found from the original image.

Figure 8 thus presents the mapping of image 11 onto image 8. This mapping is issued from the several motion fields computed between these two moments. Contours of the original image 8 are overlayed on it.

A Deformable Structure: The Heart

During the beating of the heart we observe three phenomena: the blood fluid dynamic, the cardiac wall deformations and the interaction between the blood and the cavities. Our approach is dedicated to the numerical estimation of the wall deformations because we use at once a gradient information and a curvature information which are pertinent only near contour points.

2D Case. The 2D sequence of the heart motion was acquired by ventriculography with catheterization. This technique consists in introducing an X-ray contrast agent inside the left ventricle. The contrast agent allows to evaluate the motion of the left ventricle wall. Figure 1 shows one image produced by this device. The size of the images is 256 × 256. The cardiac cycle is sampled with 38 intermediate phases, which means that the time resolution is better than for the knee sequence. Thus, the maximal displacement between two successive frames is about six pixels. For our processing we use 46 × 46 meshes. Figure 1 also shows our mesh adapted to the contours of the image.

Figure 9 shows two motion fields computed on these data. One is the motion field computed between two successive frames. Mapping results are next presented in Fig. 10. We present the mapping of images 3, 6 and 9 onto image 2. They are presented twice without and then with contours of the original image 2.

We make use of analytical modal analysis as an instantaneous method for analysis and smoothing of the computed displacement field. This technique consists to decompose our output (displacement vectors field) in the analytical basis of the free vibration modes (for more details see [START_REF] Nastar | Spatio-temporal analysis of nonrigid motion from 4D data[END_REF]). The use of modal analysis is interesting first to characterize a motion field with a few number of parameters, next to quantitatively compare several motion fields.

The Fig. 11 shows an example of modal analysis. The original motion field corresponds to the overall contraction of the heart. We show the first two modes and the motion obtained by addition of the fifty more representative modes. With only 50 modes (out of 2116) the mean error between the original field and the approximated field is about one pixel. This means that the motion information may be compressed with a factor 44 without loss of accuracy.

3D case. We make use of a set of 4D (3D images plus time) data of a beating canine heart (see an alternative use of the same data in [START_REF] Nastar | Spatio-temporal analysis of nonrigid motion from 4D data[END_REF]).

The images, obtained by fast cine-CT, come from the Dynamic Spatial Reconstructor (DSR) machine. It provides real-time data of the overall cardiac cycle without using synchronization, which is not the case until now for MRI. The size of the 3D images is 98 × 100 × 110. We did experiments on the systole (heart contraction) which is divided into seven moments. We use 3D adaptive meshes of size 25 × 25 × 25. Concerning the computational time, it takes about one minute to compute a motion field between two successive images with a workstation Dec Ultrix. This time only refers to the resolution of the Euler-Lagrange equation, and not to the preliminar computations, e.g., differential features and mesh adaptation. Indeed, since (i) our algorithm is O(N ) with N = size of the stiffness matrix and (ii) the size of the stiffness matrix is reduced in that case by a factor 60, we can estimate that it would take more than one hour to compute the motion on each voxel with the same technique. Others techniques like octree-splines could certainly be used to reduce further the computing time [START_REF] Szeliski | Matching 3-D anatomical surfaces with non-rigid deformations using octree-splines[END_REF].

First we show in Figs. 12 and 13 our motion field corresponding to the overall heart contraction. For display purposes, these views are presented with two surfaces of iso-intensity, one is the endocardium (internal wall) at the beginning-of-systole, the other is the endocardium at the end-of-systole. For a description of the algorithm used to extract the iso-intensity surfaces, one can refer to [START_REF] Thirion | The 3D marching lines algorithm and its applications to crest lines extraction[END_REF].

For checking the accuracy of this motion field, we again present some image mapping results. We first show the endocardium contours at the beginning-ofsystole overlayed on the end-of-systole image. It helps to evaluate the amplitude of the deformation during the systole. Next we show the mapping of the end-ofsystole image into the beginning-of-systole image. We also overlayed on it the original endocardium contours at the beginning-of-systole. Since these are 3D images, we only show a few slices (Fig. 14).

We finally present some results concerning the modal analysis of our motion fields. In that case, it is too long to compute all the modes and next to classify them with respect to their relevance. However, it is possible to approximate the original field just by considering the lowest frequency modes. Table 2 presents the mean error (in voxel) between the original field and a few modal approximations. Thus, with only 343 modes (out of 15625), the mean error is already less than one half voxel.

The last result illustrates one potential capability of modal analysis. Because this technique allows to describe a field with a few number of parameters, it is possible to compare efficiently different kinds of deformation. In Fig. 15, we show three spectra which represent the modes representativity as a function of the modes frequency, for each direction. The mode representativity simply is obtained by computing the scalar product between the mode and the original field. The spectra of two different deformations are printed in the same graphic.

Conclusion

We have described a method for computing dense nonrigid motion fields in time sequences of 2D or 3D images. It consists of minimizing an energy which takes into account differential similarities for contour points like gray level gradient and curvature and imposes a regularity constraint for non-salient points. In order to reduce the numerical complexity, we have introduced adaptive meshes. Since the number of nodes is less than the original number of pixels (or voxels), the computational time is drastically reduced. However, because the nodes are concentrated in the relevant areas of the image, the result accuracy is preserved. Some preliminary results concerning a quantitative analysis of the deformations are performed via a straightforward analytical modal analysis [START_REF] Nastar | Spatio-temporal analysis of nonrigid motion from 4D data[END_REF][START_REF] Benayoun | Dense non-rigid motion estimation in sequences of 3D images using differential constraints[END_REF]. Directions for further work include improving our mesh model. We intend to use unstructured meshes [START_REF] Vasilescu | Adaptive meshes and shells[END_REF]. This technique allows to reduce strongly the number of nodes and also to speed up the computation of motion. We also consider introducing other types of regularity constraints, for example in order to preserve some "natural" discontinuities of the motion, like those which occur at the boundary between healthy and ill tissues.

However, we think the most promising direction concerns the introduction of physical information about the motion. It can be very helpful in order to improve the accuracy of the results, in particular concerning the motion paths of individual points. Up to now we have considered only geometric features to recover the motion. By physical information, we mean for example the data provided by Tagged MRI [START_REF] Gutman | Tag and contour detection in tagged MR images of the left ventricle[END_REF][START_REF] Azhari | A noninvasive comparative study of myocardial strains in ischemic canine hearts using tagged mri in 3-D[END_REF][START_REF] Park | Volumetric deformable models with parameter functions: a new approach to the 3D motion analysis of the LV from MRI-SPAMM[END_REF]. This technique directly gives the true displacement of few 3D points during the heart systole. Because our method is really volumetric and tags are located everywhere, including such information would be quite straightforward. It would help to better recover the non-tangential part of the displacement, like the twisting component in the heart contraction. features from the partial derivatives of the image function [START_REF] Monga | Using partial derivatives of 3D images to extract typical surface features[END_REF][START_REF] Monga | Using partial derivatives of 3D images to extract typical surface features[END_REF]. We just recall here the main results.

Partial derivatives of the images are computed by application of a convolution filter. We used the Canny-Deriche filter which is optimal for detecting step edge contour model [START_REF] Monga | Recursive filtering and edge closing: two primary tools for 3D edge detection[END_REF]. For 1D input data, we recall here the smoothing operator corresponding to this filter:

F(x) = c 0 (1 + α|x|)e -α|x|
where α controls the width of the filter and c 0 is a normalization coefficient depending on α.

This filter has been extended to provide first, second and third order partial derivatives of the 3D images [START_REF] Malandain | Filtrage, topologie et mise en correspondance d'images médicales multidimensionnelles[END_REF].

A.2. Binary or Floating Features

For our experiments we used two categories of energies. Concerning the gradient dependent energy, we have considered two kinds of features:

• the gray level gradient norm • the "contourness" criterion which is defined as the locus where the gradient norm is maximal in the gradient direction

The choice of feature may give appreciably different motion fields. For example, a contour point may have low gradient norm and consequently will be better tracked with the contourness criterion. For the curvature dependent energy we also have tested two features:

• the principal(s) curvature(s) • the "cornerness" criterion, like the curvature extremality.

We now recall the mathematical formula which give these features. We focus on curvature(s), formula concerning contourness are given in [START_REF] Monga | Recursive filtering and edge closing: two primary tools for 3D edge detection[END_REF].

A.3. 2D Features

We suppose that a curve C is locally defined as an isointensity line. Consider C(s) = (x(s), y(s)) a normal parameterization of this curve:

I (C(s)) = I 0 , ∀s
where I 0 is the iso-value. If we differentiate twice this expression with respect to s, we obtain an expression of the curvature K of the curve depending only on the partial derivatives of the image function I :

K (C(s)) = -t T Ht g (A1)
where t is the unit curve tangential vector, g is the curve normal vector and H is the Hessian matrix of I .

By differentiating again this expression with respect to s, we have an expression of the curvature extremality E (the curvature extremums corresponding to the loci where this expression cancels):

E(C(s)) = 3 (t T Ht)(g T Ht) -g 2 t T ∂ H(t)
where ∂ H is the vector operator:

∂ H (t) =      t T ∂H ∂ x t t T ∂H ∂ y t     

A.4. 3D Features

We consider here an iso-surface S(u, v) = (x(u, v), y(u, v), z(u, v)):

I (S(u, v)) = I 0
where I 0 is the iso-value. On a point P ∈ S, we consider an orthonormal basis (h, f ) of the tangent plane T P . By a similar computation than in the 2D case, we obtain the expression of the normal curvature K d where d ∈ T P :

K d = -d T Hd g
The principal curvatures (K i ) i=1,2 correspond to the highest and to the lowest value of K d :

K i = h T Hh + f T Hf 2 g ± (h T Hh -f T Hf ) 2 + 4(h T Hf ) 2 2 g
We simultaneously obtain the corresponding principal curvature directions:

e i = h + λ i f , i = 1, 2 (A2)
with

λ i = g K i -h T Hh f T Hh
This also gives the mean curvature (K m ) and the Gaussian curvature (K g ):

K m = h T Hh + f T Hf 2 g K g = (h T Hh)(f T Hf ) -(h T Hf ) 2 g 2 .
By differentiating a principal curvature in the associated principal curvature direction, we finally obtain an expression of the curvature extremality:

E i = ∇ K T i e i = 3 e T i He i (g T He i ) g 2 e T i ∂ H(e i ) + 2α i g e T j He i with α i depending on the cross-extremality ∇ K T i e j :

α i = 1 K j -K i ∇ K T i e j

Figure 1 .

 1 Figure 1. 2D adaptive mesh on the ventriculography data.

Figure 2 .

 2 Figure 2. 2D adaptive mesh on the knee data.

Figure 3 .

 3 Figure 3. 3D adaptive mesh on the dog heart sequence (we cut non-relevant nodes for clarity).

Figure 4 .

 4 Figure 4. High pass filter (v min = 0, v max = 20).

Figure 5 .

 5 Figure5. This figure shows the evolution of a 2D curve at two successive times. The high curvature point P is not attracted by the closest contour point Q but, thank to the weighted energies defined in the text, is attracted by the closest high curvature point R.

Figure 6 .

 6 Figure 6. This figure shows the evolution of a 2D curve at two successive times. The curvature extremum point P is first "attracted" by a similar curvature point Q. Next, by relaxing the curvature constraint it matches the closest contour point R.

Figure 7 .

 7 Figure7. The top images are the gray level data of images 8 and 9 in the knee sequence. On bottom left is the mesh adapted to image 8. On bottom right is the motion field from image 8 to image 9 computed on each node of the adaptive mesh.

Figure 8 .

 8 Figure 8. The image on the left is the gray level data of image 11 (bottom). The contours derived from image 8 are overlayed on it to show the motion amplitude. The image on the right represents the mapping after resampling of image 11 onto image 8. The contours found from the original image 8 are again overlayed as a check for accuracy in the mapping estimation.

Figure 9 .

 9 Figure 9. Motion field from image 5 to image 6 computed on each node of our adaptive mesh.

Figure 10 .

 10 Figure 10. Third line: mapping after resampling of images 3, 6 and 9 onto image 2. Fourth line: same results but contours of the original image 2 are overlayed.

Figure 11 .

 11 Figure 11. Modal approximation on 2D heart data. From up to bottom and left to right: the original motion field corresponding to the heart contraction, the two more representative modes and the field approximated by addition of the fifty more representative modes.

Figure 12 .

 12 Figure 12. Profile view of the heart systole. The endocardium iso-intensity surface is grey at the beginning-of-systole and red at the end-ofsystole.

Figure 13 .

 13 Figure 13. Back view (view from apex) of the heart systole.

Figure 14 .

 14 Figure 14. On left, a few slices of the original end-of-systole image. Endocardium contours at the beginning of systole are overlayed on them. On right, the corresponding slices of the mapping of the endof-systole image into the beginning-of-systole image. Endocardium contours at the beginning of systole again are overlayed as a check of accuracy for the mapping. The mapping preserves the basic shape of the endocardium.

Figure 15 .

 15 Figure 15. Modal Directional spectra: Modes representativity versus modes frequency in each direction. Two deformations are presented : a low deformation from beginning-of-systole to intermediate step during the systole (black), a high deformation corresponding to the overall systole (grey).

Table 1 .

 1 Synthetic example: Statistics on the distance between the displacement obtained with our algorithm and the exact displacement given by the affine transformation.

			High gradient High curvature
	Distance	All nodes	nodes	nodes
	(Voxel)	(1331 Nodes) (427 Nodes)	(147 Nodes)
	Minimum	0.0321	0.0321	0.0619
	Maximum	1.4400	1.4400	0.8510
	Mean	0.4720	0.4190	0.4380
	Standard deviation	0.2380	0.1780	0.1540

Table 2 .

 2 Mean Distance Error between original motion field and some modal approximations by the frequency modes (notice that all number of nodes are power of 3 because we consider all combinations of lowest frequency modes in each direction).

	Lowest frequency	Mean approximation
	modes	error (Voxel)
	27	0.85
	64	0.75
	196	0.50
	343	0.48
	1000	0.38
	15625	0.00

Acknowledgments

Special thanks are due to Chahab Nastar for his contribution concerning modal analysis and to Professor Mike Brady for a substantial help in the writing of this paper. We wish to thank also Isaac Cohen, Alexis Gourdon and Laurent Cohen for helpful discussions. Thanks are also due to Dr. Richard Robb at Biomedical Imaging Resource, Mayo Foundation/Clinic, for providing us the fast-CT heart data. This work was partially supported by Digital Equipment Corporation. Note 1. Thanks to Dr. Boisgard from the hospital Montpied at Clermont-Ferrand (France).

Appendix

A.1. Partial Derivatives on Digital Images

Our tracking method is based on the comparison of differential features computed on the images. We already have shown in a previous study how to calculate these


Nicholas Ayache
email: nicholas.ayache@sophia.inria.fr

-Stéphane Cotin

-Hervé Delingette

Jean-Marie Clement

-Jacques Marescaux

Simulation of Endoscopic Surgery

Introduction

Medical images contain a lot of information on the internal geometry and on the functions of the human body. The introduction of local storage capacities in each production site, and the introduction of telematic links between production sites and surgical departments with high bandwidth networks have open new possibilities in the preparation, control and validation of therapy. This is particularly true for video-surgery, where the use of video-images during the intervention makes the use of medical and virtual images during training or rehearsal sessions quite natural.

The precise geometrical and physical information given by medical images can be used to create a model of an average virtual patient, in order to train the surgeon on different types of surgical procedures. Also, pre-operative images of a specific patient can be used to create virtual « copy » of a particular case, in order to rehearse a difficult surgical procedure. In the future, we believe that surgery simulators should become as common for surgeons as flight simulators for pilots.

In this article, we describe the virtual environment we have created in order to develop a simulator of liver surgery including visual and haptic feedback: the surgeon can interact in real time with a image of the virtual organ, and both look and feel the deformations. This work was done in the Epidaure project at INRIA Sophia-Antipolis (France) in collaboration with IRCAD in Strasbourg (France). Epidaure is a research group specialized in medical images analysis, and IRCAD is a multidisciplinary research institute for cancers of the digestive tract, involving among others clinicians, surgeons and scientists.

Laparoscopic Surgery

Laparoscopic surgery belongs to a set of techniques called Minimally Invasive Therapy, where the intervention in the abdomen is made through small incisions with elongated instruments. Images of the operation field are acquired thanks to an optical fiber also introduced through a small incision and linked to a mini-camera, and are displayed on a TV-screen placed above the patient (cf. figure 1).

This approach has many advantages for the patient: small incisions, reduced operation traumatises, and post-operative morbidity. Also the hospital period of stay can be drastically reduced. On the other hand, this technique requires new capabilities from the surgeon. First, degrees of freedom of the instruments are reduced, therefore introduction in the abdomen has to go through a fixed point (troccard). Second, the vision of the operation field is distorted by wide angle lenses, and does not correspond to an observation from the surgeon's point of view. In addition, the location of the TV-screen is not in the direction of the instruments. This requires a particular hand-eye coordination.

Simulation of hepatic surgery

Therefore a surgeon must practice thoroughly this new technique with the new sensations attached to it. Training on cadavers or animals raises a number of practical and ethical problems. Mechanical devices are extremely limited and do not offer a realistic simulation. The advantages of a computer simulator are numerous. Different organs and pathologies can be represented within a common framework, in order to train the surgeon to a large number of situations. Also, they can be adapted to the case of a specific patient. A number of groups (TIMC -Grenoble University Hospital -LIFL -University of Lille -France, KISMET -Karlsruhe -Germany, HT -Washington and MERL -Boston -USA) regrouping computer scientists and physicians are working towards the development of surgery simulators. The collaboration between INRIA and IRCAD has at two goals: first the creation of an anatomical model of the liver to teach anatomy; second, the development of a surgery simulator to teach medical gestures.

Geometrical modelling

In order to create a very realistic model of the liver, we used the images of the Visible Man, distributed by the National Library of Medicine (USA). These images are high resolution (1 mm) colour pictures of cryogenic cross-sections of a cadaver. We used the digital image processing tools of the Epidaure group to extract the external boundary of the liver, and also the various vessels (veina cava, portal vein, hepatic artery, billiary duct). Then a volumetric representation is created (). Finally, an isosurface extraction [START_REF] Thirion | Computing the differential characteristics of iso-intensity surfaces[END_REF] provides the mesh of a 3-D surface corresponding to each structure (cf. figure 3). One obtains a set of 5 surface meshes, either triangulated or with a simplex mesh [2] [3]. These surfaces can be manipulated interactively through a graphical interface. The surgeon can therefore plan surgery with the accurate knowledge of vessels location with respect to the liver external boundary, which is the only visible item during the actual surgery. Thanks to this 3D reconstruction, it is possible to visualise with transparency the different vessels of the liver, and therefore to plan the impact of each medical gesture. Doing this can reduce drastically the risk of cutting a major vein or artery during surgery.

Physical modelling

The creation of an accurate model of the liver is only one step for surgery simulation. In effect, as said before, a major goal is to render the sensations linked to laparascopy, in order to better practice the required hand-eye coordination. For visual realism, we must simulate what appear on the video screen. This implies that it is not necessary to be « immersed » in a virtual world for this kind of simulation. On the other hand, the tactile feed-back is an essential complement to the visual feed-back, and we believe that an efficient simulator must be equipped with force feedback.

Organs are deformable, and their deformation laws are in general complex and unknown. A number of approaches have been proposed to model deformable objects in graphics [START_REF] Terzopoulos | Elastically Deformable Models[END_REF][START_REF] Terzopoulos | Modeling Inelastic Deformation : Viscoelasticity, plasticity, fracture[END_REF][START_REF] Platt | Constraint Methods for Flexible Models[END_REF][START_REF] Bro-Nielsen | Modelling elasticity in solids using active cubes -application to simulated operations[END_REF], or in a domain closer to physiology [START_REF] Waters | A Physical Model of Facial Tissue and Muscle Articulation[END_REF] [START_REF] Kuehnapfel | CAD-Based Graphical Computer Simulation in endoscopic surgery[END_REF]. We have chosen a relatively simple model, but with a lot of advantages including speed. Our objective is to develop a simulator which does not remain a laboratory prototype, but could be industrialised for surgeons benefit. To achieve this goal, the chosen material must be standard and of reasonable cost. Our chosen model should allow:

• real time deformations (25Hz for display, between 300Hz and 500Hz for force feedback).

• volumetric modeling of organs, allowing realistic deformations and simulation of cutting, suturing, etc. • modelling anisotropic linear elasticity properties of liver tissues and vessels (although their is a lack of biomechanical data available on this topic).

Force feedback

Recent work has shown that the notion of presence in a virtual environment is not only dependant on visual realism but greatly enhanced by other perception means [START_REF] Gibson | Simulating Arthroscopic knee suergery using volumetric object representation, realtime volume rendering and haptic feedback[END_REF][START_REF] Barfield | Factors Affecting Presence and Performance in Virtual Environments[END_REF]. In particular, the sensation of forces and textures is very important in medical applications [START_REF] Sukthankar | Force Feedback Issues in Minimally Invasive Surgery[END_REF][START_REF] Jackson | Force Feedback and Medical Simulation[END_REF]. Force feedback remains essential to simulate realistically the technique. This is the only way to restore the « physical » contact between the operator and virtual objects.

We chose to use the force feedback system called Laparoscopic Impulse Engine™. This system has many attractive properties including 5 degrees of freedom (3 of them being motor controlled), an accuracy in displacement of about 0.025mm and maximum intensity force in the direction of the trocard of 8.9N. This system is integrated in an architecture with 2 workstations linked with ethernet (cf. figure 6 and figure 11) .

To integrate such a system with our simulator, it is necessary to determine which forces are implied. To respect the constraints of speed and realism, linear elasticity theory appears rather attractive. This theory allows a better physical modelling than mass-spring systems, in particular because the reacting force is linked to the object deformation and to its elastic properties by the continuum mechanics relations. Linear elasticity has often been used [START_REF] Jackson | Force Feedback and Medical Simulation[END_REF][START_REF] Thomas | Three Dimensional Finite Element Analysis of Elastic continua for Tactile Sensing[END_REF][START_REF] Song | Tissue Cutting in Virtual Environment[END_REF] as a good approximation -for small displacements -of a deformable body behaviour, but the resolution of the attached equations requires finite element techniques which are quite costly in terms of CPU. To reach real time performances, we have introduced specific pre-computations of elementary deformations.

Force feedback system Control unit (PC) Workstation with 3D graphics capabilities

Figure 6: the actual hepatic surgery simulator prototype includes a workstation with 3D graphics capabilities, a control unit (PC Pentium) and a force feedback system « Laparoscopic Impulse Engine ».

3D linear elasticity

Let us recall briefly the bases of linear elasticity. Let Ω 0 be the configuration of an elastic body before deformation. Under the action of a field of volume and surface forces, the elastic body deforms itself and takes a new configuration Ω. The problem is then to determine the displacement field u which associates to the final position P 0 of each particle of the body before deformation, its position P in the final configuration.

In order to find a numerical solution to this problem, we used finite elements through a variational formulation. The problem finally reduces to the resolution of a linear system

[K] [U] = [F]

where [K] is the stiffness matrix describing the elastic properties of the solid, [U] corresponds to the unknown displacement field and [F] represents the external forces. We chose a decomposition of the object volume in a set of tetrahedra, because we already have a triangulation of the outer boundary of the liver. In addition , the geometry of most organs being rather complex, tetrahedra tend to be a good trade-off for such a decomposition.

In addition to the algorithmic complexity, another problem linked with the of finite elements techniques is due to boundary conditions. In a classical approach, one considers that the deformation is induced by a set of forces applied to the surface of the model and combined with a set of boundary conditions (usually a set of nodes with a null displacement). The force feedback system does not allow the measurement of forces generated by the operator, but only the measurement of the displacement corresponding to the mobile part of the force feedback system. Therefore, we must solve a problem in different format, where our boundary conditions include the displacements of the organ nodes which are in contact with the instrument tip. From this, we have to deduce the global deformation of the system, and the reaction force finally sent back to the operator through the system. (cf. figure 11).

Real time computation of the deformation and forces

Because of high frequency computations required by force feed-back systems, we have developed a method allowing the computation of both the global model deformation and of the reaction force at a frequency varying between 300 and 500Hz.

Our algorithm is based on the principles of linearity and superposition of elasticity theory. The superposition principle can be summarised as follows: if we apply to the solid a set of forces creating a displacement field u, then this same displacement field can be obtained by applying separately each force, and by adding up their effects. The linearity principle tells us that if a force f, applied to the surface of an elastic body creates a displacement u, then a force 2f will create a displacement 2u.

Consequently, any deformation of our elastic model can be obtained by a linear combination of a set of « elementary » displacements applied to each surface node. An adaptation is required in our case, because displacements are imposed, not forces [START_REF] Cotin | Real-time Volumetric Deformable Models for Surgery Simulation[END_REF][START_REF] Cotin | Geometric and Physical Representations for a Simulator of Hepatic Surgery[END_REF]. Our method allows the real-time computation of the deformation of a complex organ like the liver under the action of a virtual instrument, along with the reaction force that this organ transmits to the surgical tool (cf. figure 7, figure 8 and figure 9).

Perspectives

It is not sufficient to allow the surgeon to manipulate and deform the anatomical model of the liver, we must also ensure that the model behaviour corresponds to anatomical reality. The elastic properties, boundary conditions must be determined for the liver to allow a more realistic behaviour. Visual realism can also be enhanced. We propose the following future improvements:

• nature of hepatic tissues: we attach a elasticity coefficient which depends on the grey-level value of the medical image of the liver. Incidentally, we could take into account the vascular rigidity which plays an important role in the general behaviour of the liver.

• interactions with other abdominal organs: the liver is not isolated in the abdomen, but it is surrounded by other organs like the gall bladder. Consequently, it is necessary to model these interactions between nearby organs, which could be done by an additional set of boundary conditions.

• texture: a realistic texture, extracted from real video-surgery images could be mapped on the surface of virtual organs to enhance the visual realism.

Conclusion

This work must be considered as a first stage towards the building of a practical surgical simulator. We have demonstrated that real time deformation of an elastic model was possible, by combining a physical model with an adequate geometrical representation. We believe that it is essential to have a volumetric modelling of organs, in order to produce realistic deformations, and allow, in the future, the simulation of a realistic hepatectomy. Developing an original and efficient method for the computation of deformations based on linear elasticity, we have been able to simulate real-time interactions with a virtual organ and to compute the contact forces with the surgical instrument. These forces have a physical foundation as they depend on the imposed displacement and are based on linear elasticity. The prototype we have developed allowed for the validation of the algorithms, and a better understanding of the software integration constraints. A new version of this simulator should migrate from INRIA to IRCAD in the near future for further evaluation steps.

Figure 1 :

 1 Figure 1: Minimally Invasive Surgery principle. The surgeon can see the operative field using a micro video camera. Consequently, the manipulation of the surgical tools is more complex than conventional open surgery. Origin : Forschungszentrum Karlruhe

Figure 2 :

 2 Figure 2:. liver model built from a set of 2D segmented medical images.

Figure 3 :

 3 Figure 3: result of the vascular tree segmentation. The high quality of the original images allows a high level of detail in the reconstructed model.

Figure 4 :

 4 Figure 4: surface mesh of the liver. A minimisation of the number of vertices -without a lack of accuracy -is possible with the use of an adaptive simplex mesh.

Figure 5 :

 5 Figure 5: anatomical liver model. The visualisation of the vascular tree is possible with the use of transparency properties. It is possible to manipulate the 3D mesh in real-time.

Figure 7 :

 7 Figure 7: the liver model in its rest shape. The square landmarks correspond to the boundary conditions.

Figure 8 :

 8 Figure 8: first, the surgical tool comes in contact with the mesh surface. Then, its following motion is converted into a displacement of a set of nodes on the surface.

Figure 9 :

 9 Figure 9: the motion of the surface nodes provides a deformation, computed in real-time. Moreover, a force associated to the deformation is also computed and sent to the force feedback system.

Figure 10 :

 10 Figure 10: image sequence: real-time interactions with a complex liver model are described. The gall-bladder is also represented. The surgeon feels a force which is correlated with the deformation amplitude and with the geometry of the mesh.

Figure 11 :

 11 Figure 11: the actual hepatic surgery simulator prototype at INRIA. At the fore-front we can see the force feedback system. The screen shows a scene viewed from the endoscope view point.