Nicholas Ayache

Keywords: Medical computer vision, virtual reality and robotics volume image processing, medical images, 3D vision, virtual reality, medical robotics, research trends

The automated analysis of 3D medical images can improve both diagnosis and therapy significantly. This automation raises a number of new fascinating research problems in the fields of computer vision, graphics and robotics. In this paper, I propose a list of such problems after a review of the current major 3D imaging modalities, and a description of the related medical needs. I then present some of the past and current work done in our research group EPIDAURE* at INRIA, on the following topics: segmentation of 3D images; 3D shape modelling; 3D rigid and nonrigid registration; 3D motion analysis; and 3D simulation of therapy. Most topics are discussed in a synthetic manner, and illustrated by results. Rigid matching is treated more thoroughly as an illustration of a transfer from computer vision towards 3D image processing. The later topics are illustrated by preliminary results, and a number of promising research tracks are suggested.

INTRODUCTION New images

Three-dimensional (3D) images are becoming very popular in the medical field'-'. A modern hospital each year commonly produces tens of thousands of volumetric images. They come from different modalities like magnetic resonance imagery (MRI), computed tomography imagery (CTI, also called scanner imagery), nuclear medicine imagery (NMI) or ultrasound imagery (USI).

These images share the particularity of describing the physical or chimical properties at each point of a *EPIDAURE is a very nice location in Greece which used to be the sanctuary of ancient medicine. Today, for computer scientists, it is also a recursive acronym (in French): Epidaure, Project Image, Diagnostic AUtomatique, et RobotiquE.

INRIA-EPIDAURE

Project, BP 93-06902 Sophia-Antipolis, France Paper received: 8 March 1994; revised paper received: 7 February I995 studied volume. This information is stored in a discrete 3D matrix Z(i, j, k) of voxels (volume elements), called a 3D image because of the analogy with digital 2D images Z(i, j) stored as 2D matrices of pixels (picture elements).

Actually, the most simple exploitation of 3D images consists in a direct visualization of the successive 2D images (also called cross-sections) obtained by setting one of the three coordinates (i, j, k) to a constant. In fact, applying advanced 3D image processing and graphics allows much more vivid and realistic displays, as shown in CT1 and MRI provide complementary anatomical information.

CT1 gives very contrasted images of the bones, while MRI gives excellent descriptions of most organs and soft tissues, as well as a reasonably good description of the bones.

NM1 measures locally the density of some injected radioactive molecules.

It therefore provides functional information on the metabolism of the regions studied (e.g. the density of glucose used per unit of time and volume in the brain). The main sources of NM1 are single photon emission computed tomography (SPECT), with positron emission tomography (PET) still being a research tool.

Finally, US1 measures locally the variation of the acoustic impedance (a function of the speed of ultrasound propagation).

Although most commonly used in two dimensions, ultrasound images can also be acquired in 3D", ". Ultrasound images are easily acquired at a very fast rate, giving both anatomical information (boundaries of anatomical structures usually correspond to variations of the acoustic impedance), and dynamic information (e.g. the heart motion). Dynamic information can also be obtained by gated (i.e. synchronized) MRI or NMI.

Although US1 and MRI are recognized as noninvasive (non-dangerous) techniques, this is not the case for CT1 and NMI, which involve radioactive materials. The cost of equipment varies a lot from one modality to another: currently the cheapest is US1 (typically, less than 100 k$), while the most expensive is MRI (typically several M$).

New 3D imagery devices are emerging, like angiographic MRI, which describes the anatomy of the vascular system, magneto-encephalography equipment, which measures magnetic field variations, or functional MRI, which provides metabolism information (like NMI), but with a non-invasive method.

These new modalities are still at the research stage, but may lead to major 3D imaging devices in the near future.

New markets

The market of the production of medical images was evaluated at being 8 billion (US dollars) in 199112"3, and shows approximately an increase of 10% per year. Among these figures, MRI represents currently a market of 1 billion, with a strong increase of approximately 20% per year14.

Besides the production itself, 3D image processing is the most recent market. Almost non-existent a few years ago, it was evaluated at $350 million in 1992, with a planned evolution between 2&40% per year during the next five years.

New medical needs

Exploiting 3D images in their new format (a 3D matrix of numbers) is usually a very awkward task. For instance, it is now quite easy to acquire MRI images of the head with a resolution of about a millimeter in each of the three directions.

Such an image can have 2563 voxels, which represents about 17 Mbyte of data. High resolution 3D images of the heart can also be acquired during a time sequence, which represent spatiotemporal data in four dimensions.

In both cases, displaying 2D cross-sections one at a time is no longer sufficient to establish a reliable diagnosis or prepare a critical therapy.

Moreover, the complexity of the anatomical structures can make their identification difficult, and the use of multimodal complementary images requires accurate spatial registration.

Long-term study of a patient's evolution also requires accurate spatio-temporal registration.

We identified the automation of the following tasks as being crucial needs for improvment in diagnosis and therapy:

1. Interactive visualization must really be 3D, with dynamic animation capabilities.

The result could be seen as a flight simulation within the anatomical structures of the human body. A recent review of the state-of-the-art can be found in Stytz et aL2'. 2. Quantifi:cation of shapes, textures and motion must provide the physician with a reduced set of parameters useful to establish his diagnosis, study temporal evolution, and make inter-patient comparisons. This must be true for the analysis of static and dynamic images. but also because usual approximations like polyhedral models or typical assumptions like rigidity rarely apply to medical objects. This opens a large range of new problems sometimes more complex than their counterparts in 2D image analysis. On the other hand, specific properties of 3D medical imagery can be exploited very fruitfully. For instance, contrary to video images of a 3D scene, geometric measurements are not projective but Euclidean measurements. Three-dimensional coordinates of structures are readily available! Moreover, one can usually exploit the intrinsic value of intensity, which is generally related in a simple way to the physical or physiological properties of the considered region; this is almost never the case with video images, where intensity varies with illumination, point of view, surface orientation, etc.). Also, a priori knowledge is high, in the sense that physicians usually have protocols (unfortunately, depending on the image modality) to acquire images of a given part of the body, and different patients tend to have similar structures at similar locations! Finally, having a dense set of 30 data provides a better local regularization when computing local differential properties, as we shall see later.

New computer vision, graphics and robotics issues

Having listed the medical needs and the new image properties, I now set a list of computer vision, graphics and robotics issues which I believe are central problems: It is necessary to describe non-polyhedral 3D shapes with a reduced number of intrinsic features. This involves mainly computational and differential geometry. 30 matching of 3D shapes: once segmented and modelled, new algorithms must be designed to reliably and accurately match such representations together, both in rigid and nonrigid cases. This is necessary to solve registration and identification needs.

We have found that a set of generic tools was quite effective in solving completely a number of specific segmentation problems, both in static and dynamic images. These include 3D edge extraction, 3D digital topology and mathematical morphology operators, 2D and 3D active contours, and 3D texture analysis tools. We found that, combined together, these tools could solve the problem of segmenting major anatomical structures in CT1 and MRI, and also the tracking of 2D and 3D structures of the beating heart in dynamic US1 and NM125. 30 motion analysis: this requires the development of new tools to process sequences of 3D images (i.e. 4D images!), in order to track and describe rigid and nonrigid motion of 3D anatomical structures. This is necessary for the quantification of motion needs. Dynamic physical modes of anatomical structures should be developed to provide realistic simulation of interaction with 3D images. This is required for the planning and simulation of therapy. Physical models can also help solving the previous 3D motion analysis problems. Geometric reasoning is required to help therapeutic planning, in particular to determine trajectories of beam sources in radiotherapy, and the succession of accurate medical gestures in surgery. Virtual reality enviranments should be developed to provide realistic interactive visualization and to help planning and simulation. Dedicated medical robots, possibly passive or semiactive, equipped with specific sensors (force sensing, optical or ultrasound positioning, etc.) must be developed for the automatic control of therapy.

3D edges

As previously mentioned, a major advantage of 3D images comes from the fact that intensity is usually a simple function of the structures studied (contrary to video images of a 3D scene). To segment CT1 or MRI images, assuming constant intensity plus additive noise within a given anatomical structure is often a reasonable assumption in regions with low texture. Intensity in NM1 images is also simply related to the physiological function studied, but the resolution is often lower with a higher noise level. The most difficult images to segment are probably USI, where intensity already measures a function of the derivative of the acoustic impedance of the tissues, but with a strong multiplicative (in frequency domain) noise producing a very typical 'speckle' texture almost everywhere.

For regions with low texture, Monga26'27 showed that the 3D generalization of the Deriche-Canny edge operato?

was quite efficient to extract edges. The superiority of 3D filtering the volumetric data instead of successively tiltering 2D cross-sections was clearly demonstrated, as can be seen in

SEGMENTATION OF 3D IMAGES

Because

3D filtering is computationally more intensive, the use of separable recursive filters is crucial for the sake of computational complexity. Malandain implemented a version of the 3D edge detector which has been distributed to several places.

Digital topology and mathematical morphology

Thresholding and/or edge detection must generally be followed by some 3D mathematical morphology

Deformable surfaces

This approach consists in generalizing in three dimensions the active contours introduced by Terzopoulos and coworkers36737. A variational method has been developed by Laurent Cohen, and then by Isaac Cohen, to minimize an energy which accounts for the elasticity properties of the surface (internal energy) and the proximity of detected edgels (external energy). The mathematical framework is one of the finite elements, which provides a sound and stable analytical description of the surface. This approach provides a good segmentation of surfaces with a simple topology3smm4'. We show in Figure 6 an example of the segmentation of a heart ventricle in a 3D image. After the segmentation of such a 3D image, deformable surfaces appear to be sufficiently robust to track a deformable structure in a time sequences of 3D images4'. Connected work can be found elsewhere42ms2.

Sonar space filtering

In some cases, it might be interesting to perform the segmentation before 3D image reconstruction, i.e. from the raw data obtained by the original sensor. This was clearly demonstrated for USI, with a method called sonar space filtering.

In USI, Herlin and Vojak showed that the polar geometry of US1 requires a special type of filtering which computes edges at the local resolution of the raw data. This is important, because once the image is rectified from polar to Cartesian coordinates,

Geometric tomography

In the same spirit as above, but for CTI, Thirion showed that it was possible to extract the boundary of contrasted objects directly from the edges of the sinograms (projections from which tomography is computed).

Texture-based approaches

Levy-Vthel, Berroir and Mignot58-60 implemented a system called ARTHUR, which combines texture modelling and a sophisticated discriminant analysis scheme to select a set of discriminant texture parameters from a training set of images (this approach has similarities with the work of Kubler and Gerig6'). The ARTHUR system can then be connected to a texture segmentation program called EXCALIBUR, which exploits the training phase and the notion of 'mixed classes' to provide automatic segmentation of textured images. The texture operators include 2D and 3D texture operators; some of them involving advanced fractal and multifractal parameters.

Figure 7 shows the result of such a textural segmentation in a cross-section of an MRI, at the level of the heart.

MODELLING FREE-FORM SURFACES

Once objects are segmented, they have to be modelied for further processing.

This modelling is necessary for registration (next section), which is either rigid (some This is a very challenging problem in medical image analysis, to study, for instance, the evolution of a pathology such as multiple sclerosis (MRI images), or the efficiency of a treatment against a cancer lesion (MRI or scanner images). This is also necessary to compare patients together and help therapy planning.

Extracting ridges

The approach successfully developed in the EPIDAURE project consists in extracting typical lines, called ridges or crest lines, on the surface of the objects in each volumetric image.

Our definition of ridges on a surface is the locus of points where the maximum curvature (in absolute value) is extremal in the direction of maximum curvature. This definition does not apply when principal curvatures have the same absolute value, in particular at umbilics. Anyhow, as can be seen from experimental studies (see below), these ridges convey intrinsic, meaningful and compact information on many anatomical surfaces: a intrinsic, because ridge points are invariant with respect to rigid displacements (this comes from the intrinsic properties of the curvature); 0 meaningful, because ridges tend to correspond to anatomical features (this is the case on the brain, the face and the skull, for instance, as can be seen in One should note that our definition of ridges is intrinsic to the isophote surface, without privileging any particular direction of the 3D space. One should refer to Koenderink76 for a rigorous and historical analysis of oriented ridges, i.e. ridges depending on a privileged direction in space (e.g. gravity).

Medical computer vision, virtual reality and robotics: N Ayache

Computing

curvatures requires the extraction of the differential properties of the surfaces up to the 2nd order, while ridges require 3rd order differential properties. We show in the following subsections the different approaches we investigated to achieve a reliable ridge line extraction.

Local fitting of quadrics

After a segmentation of the 3D images, using, for instance, 3D edge detectors, Monga and Sander showed that it was then possible to use a local Kalman filter around each edgel, to fit a quadric surface in the least square sense. This polynomial approximation allows the computation of differential surface properties up to second order only (curvatures).

This approach takes into account the uncertainty attached to the edge localization and the surface normal, and provides an explicit estimation of the uncertainty attached to the computation of the principal curvatures. Unfortunately, although parallel in nature, the method is computationally very expensive. Curvature information is accurately extracted, but this is not the case for ridges77.

Deriving the image intensity

Another approach is necessary to extract ridges. As mentioned earlier, it is often a reasonable assumption to say that the boundary of an anatomical surface corresponds to an isophote (or iso-intensity) surface in the volume image, especially with CTI, but also with MRI and sometimes NMI. This assumption, combined with the implicit function theorem, allows the computation of the differential properties of the anatomical surface from the derivatives of the intensity function. This method provides excellent results in high resolution volumetric images, for a computational cost which can be maintained at a reasonable limit thanks to separable recursive filtering, as shown by Monga and Benayoun63,64 and Thirion65.

Figure 8 shows the extraction of ridges on the brain surface acquired with MRI, while

Using B-spline approximations

The previously mentioned deformable surfaces had two drawbacks: high computational complexity, and loworder derivability.

The introduction of B-splines with high order polynomials, thanks to their nice separability properties, reduces these drawbacks drastically while keeping the advantages of simultaneously providing a This is very useful to track time-varying structures with rigid or nonrigid motion.

The only limitation comes from the fact that the surface topology must be known in advance. This is described in detail by Gutziec78,79.

Figure 9 shows the segmentation of a face with a deformable B-spline within an image of the head acquired with MRI, while

Marching line algorithm

In the two previous approaches (filtering and B-splines), the connection of extracted ridge points to form ridges can be achieved by a very smart algorithm invented by Thirion, and called 'the marching line algorithm'. This algorithm looks for the zero-crossings of an extremality criterion (here the maximum curvature) in an isophote surface defined by a constant intensity level I.

RIGID 3D REGISTRATION

I now describe in more detail the problem of rigid 3D registration, because we believe it is a remarkable illustration of computer vision applied to 3D medical images. Artificial landmarks could be used to simplify the point extraction and matching process. For instance, a stereotaxic frame can be rigidly attached to the patient's skull. In fact, this is not a very comfortable solution for the patient, and anyway, such a frame cannot be worn for a long period (e.g. six months!). Moreover, it can happen that some displacement occurs between the frame and the patient between the two acquisitions, or that some internal organ moves with respect to the external frame (e.g. a slight motion of the brain with respect to the skull). Finally, the accurate localization of specific points with artificial markers is usually not an obvious task.

Importance of the problem

For all these reasons, we believe that a fully automatic registration procedure relying only on detected anatomical landmarks is much more flexible and powerful than a manual procedure, or than a procedure relying on artificial landmarks.

The output of the registration procedure must be the six independent parameters of the rigid displacement undergone by the region of interest between the two acquisitions.

More precisely, one must compute the rotation and translation parameters which best match the two acquired images of this region of interest. At this point, it is important to note that because of potential occlusions, the two 3D images cannot be registered globally with a method based on the registration of the centres and axes of inertia.

Using ridges and geometric hashing

The idea is to extract the maximum number of Euclidean invariants computed exclusively on ridges, and to use them for registration.

Euclidean invariants on surface curves

We consider not only the curvature k and torsion z of the ridge lines (the parameters which characterize completely a curve up to a rigid displacement), but also the maximum curvature of the surface, ki, the angle 8 of the surface up to the 3rd order, the values of the maximum curvature value kl and direction el, as well as the surface normal N, are readily available.

But to compute the five intrinsic parameters (k, z, kl, 8, 4), we also need to compute the differential properties of the curve itself up to the 3rd order. Gukziec found that an efficient method was to approximate each ridge line (a discrete set of ridge points connected by the marching line algorithm) by a constrained B-spline" (the B-spline is constrained to have a tangent orthogonal to the surface normal). This provides a good estimation of the Frenet frame of the curve, (t, n, b), as well as the local curvature k and torsion z. It is then immediate to compute 0 = angle (n, N) and 4 = angle (t, el).

Matching algorithm

Once we know how to extract ridge points and how to compute the quantities (k, 5, k,, 8, 4) Preprocessing stage This stage can be applied offline, as soon as the first image, called the model image, is acquired.

The ridge lines are extracted from this image, and each ridge point with its Frenet frame (t, n, 6) is stored in a 5-dimensional hash table based on the five intrinsic parameters (k, z, kl, 0, 4) attached to the ridge point

The 5-dimension hash table is stored in a Idimensional array using classical but efficient hashing techniques (see Figure 13).

Recognition stage

Once the second image, called -l/F+ Buckets Real Space

Model Curves Parameter (Hash)-Space Hash Function H-Table Scene Curves

Figure 13 Preprocessing: each ridge point of the model, is stored in a hash table of dimension-5, after the computation of five independent intrinsic parameters.

At recognition time, the same ridge point, although in another position and orientation, will point towards the same cell of the hash table (courtesy of Guiziec) the scene image, is acquired, it is necessary to apply the recognition stage on-line. The algorithm considers each ridge point Si of the scene, and looks in the hash table for a model point with similar intrinsic parameters (k, r, ki, 13, 4). A difficult point is the notion of similarity, and a preliminary statistical study allowed for the definition of a statistical distance (Mahalanobis distance) based on the covariances of these parameters.

If a model point Mj is found with similar parameters in the hash table, the Frenet frames of both points are used to predict a rigid transformation T which maps both points with their Frenet frames (see Figure 14). This transformation is accumulated in a 6-dimension accumulator of rigid transformations A, by incrementing the cell whose coordinates correspond to T, and by keeping a pointer towards Fi and Si the maximum number of access. This cell contains a transformation T* with a covariance matrix C', and it is chosen as being the solution of the registration problem.

In the skull example presented, this transformation corresponds to the matching of 335 couples of ridge points. We can apply this transformation to all the ridge points of the model, and verify that they superimpose very nicely on the scene ridge points, as can be seen in Figure 16. A quantitative analysis showed that the accuracy of the registration was better than one millimeter within the volume of the head".

Matching of proteins

Pennec

found that it was possible to extend the geometric hashing ideas to the problem of finding common 3D substructures in proteins.

He found an O(n2) algorithm".

A typical result is shown in

Introduction of extremal points and extremal mesh

Thirion92T93 recently introduced a new subset of intrinsic points that he called 'extremal points', which can be identified as the subset of ridge points for which the second principal curvature k2 is also extremal in the associated direction e2. These points have the nice property of being extrema of curvature of the two principal curvatures simultaneously, and can be defined as the intersection of three iso-surfaces: a chosen isointensity surface, and the two surfaces defined by the zero-crossings of the derivatives of ki in the direction ei for i = 1,2. These derivatives can be computed everywhere but at points where the principal curvatures have identical absolute values, which again includes the umbilics.

Not only the extremal points have the property of being very stable and compact on anatomical surfaces (a few hundreds in a high resolution 3D image of the head), but some of them tend to be extremely stable from one patient to another one. They are, therefore, very good candidates as anatomical features for the Atlas Matching Problem, described in a later section.

Finally, a unifying description, called the extremal mesh, which includes ridge lines and extremal points, haas also been introduced by Thirion94,95. This description is based on the zero-crossings of the product of two extremality criterions, defined at non-umbilic points as the derivative of one of the two principal curvatures in its associated direction.

The product of these two criteria is clearly zero along ridge lines, and also at extremal points, and the definition naturally provides an intrinsic mesh on the surface, which itself proves extremely useful for its description, and also for rigid and non-rigid registration.

Matching a cross-section with a 3D image

Another interesting problem is the matching of a crosssection with a 3D image. This happens when the patient must be registered with previously acquired 3D images, during surgery, for instance.

In this case only a few (possibly a single) cross-sections are acquired with a given modality, possibly with a laser range tinder, and have to be matched with a surface extracted from a 3D image.

Relaxing the rigidity assumption is necessary in three different problems:

the atlas-matching problem; interpatient registration; and the analysis of nonrigid motion.

Atlas matching problem

Gourdon

showed that it was possible to exploit Building an 'electronic atlas' of the human body is again differential geometry constraints between the curve and a very challenging issue for the future, and several teams

Matching with a potential-based method

Computing ridges and Euclidean invariants of second and third order requires high resolution images. This is currently not possible with NMI. Therefore, in this case, and in particular to combine together images of different modality like NM1 and MRI, for instance, other methods, usually based on potentials of attractions, can be used efficiently97-'03. Contrary to the method described in the previous section, based on geometric hashing and Euclidean invariants, potentialbased methods require a preliminary estimation of the correct superposition, which might sometimes limit their robustness.

Anyhow, excellent results can be obtained, as shown in Our goal is oriented towards the automatic registration of such an atlas against the 3D image of an arbitrary patient. For doing this, we believe that ridge lines are good anatomical invariants, and can therefore serve as a sound geometrical basis to build a computerized anatomical atlas of some parts of the human body. This assumption is supported by our experimental studies, and by statistical and anthropological studies662 67.

Having built a 'generic model', the problem is then to find a matching algorithm which can register it with the 3D image of an arbitrary patient. Such an algorithm could also help in making inter-patient registrations and comparisons.

The current strategy within the EPIDAURE project, is to define a constrained network of ridge lines and extremal points92593' '06, which could deform itself to adapt to the geometry of a given patient. The spirit is

nter-patient registration problem

To register together anatomical surfaces extracted from two different patients, Feldmar introduced a quite successful approach.

He proposes to use curvature information on the two surfaces to successively find the best rigid, then globally aftine, then localiy affine transformation.

The curvature information appears to be essential to preserve the correspondance between anatomical invariants, e.g. the nose, the chin, the eyebrows . . . on a face. A typical example is shown in Figure 21. Details of this approach can be found elsewhere"2~'14.

Another approach, also quite successful, was published by Szeliski and Lavallee"'.

Analysing nonrigid motion

In time sequences of 3D images, like ,USI, MRI or NM1 images of the beating heart, it would be extremely profitable not only to track the motion of moving structures, but also to provide physicians with a reduced set of pertinent parameters which describe their nonrigid motion.

Usually, the problem is dual: first, it is necessary to find the motion of individual points; then it is necessary to approximate this motion with a reduced number of parameters.

Motion of individual points

Tracking structures with deformable models like snakes does not usually provide the motion of each individual Good results can also be obtained with a physicallybased deformable model which provides tracking and point-to-point correspondences at the same time'24,'25. Figure 23 shows the use of the deformable model of Nastar to track the mitral valve in images of the heart acquired with USI, and get point-to-point co;respondences.

Global analysis of motion

Once point-to-point correspondences have been established, it is possible to project on a reduced basis the set of displacements. This is the purpose of modal analysis, where the basis corresponds to a reduced se: af some 'qualitative' deformations126% lz7, or to the major monofrequency vibration harmonics of the elastic structure'24 (see Figure 24 for an illustrative example of the principal energy-increasing vibration modes of a cylinder). A major advantage of the latest approach is the possibility to compute analytically'*' the modes beforehand, which reduces tremendously the computational complexity. temporal Fourier coefficients, i.e. only 45 real numbers, with such an accuracy that it is possible to recreate a virtual 3D motion again almost impossible to distinguish from the real one (refer elsewhere'29-132 for a detailed presentation).

Global analysis of nonrigid motion can' also be obtained with parameterizable deformable shapes using for instance superquadrics'33m'37 or Fourier models'38. A typical example is shown in Figure 28. A necessary condition for the success of this research will be the thorough collaboration between research scientists, medical doctors and also members of related hardware and software companies.

SURGERY SIMULATION

Another condition will be the establishment of powerful electronic connections between the various departments of the hospital, as well as between hospitals.

Once these conditions are fulfilled, there is no doubt in my mind that medical image processing and medical robotics will become one of the major application fields of computer vision and robotics science, and will bring a revolution in medicine for the coming decade.

Figure 2

 2 Figure 1 3D display of two well-known anatomical structures extracted automatically from a single 3D magnetic resonance image of the head. The automatic extraction and display of these structures involves advanced 3D image processing

Figure 3

 3 Figure 3 Four different modalities of acquisition of 3D images: MRI, CTI,, NM1 and USI. Although each image is produced in three dimensions, only one cross-section is shown here. The first three images show the'head; the last is a young foetus (7 cm!)

  last one, but see18, 19,24 for spectacular examples of dedicated medical robots) by presenting the research conducted during the past five years in the research group EPIDAURE at INRIA. I present the basic lines of this research, the main results, and indicate some promising research tracks. References are made primarily to the papers published by the EPIDAURE group, although I have tried to add a significant

Figure 4

 4 Figure 4 3D edge detection is more robust and coherent than 2D edge detection.Left: result of 2D edge detection performed in successive cross-sections of the brain which were orthogonal to the plane of the picture. Result from one cross-section to the next is not as coherent as it is with 3D edge filtering (right), where none of the three principal directions is privileged (courtesy of Malandain and Monga)

Figure 5 Figure 6

 56 Figure 5 Extraction of topological singularities on the syrface of two skulls. One should note the remarkable stability of these lines, which can be used for registration (courtesy of Malandain and Bertrand)

  some information is lost, and the edge extraction results are degraded. Results obtained in time sequences of echographic images allow analysis of the motion of the mitral valve of the heart, with the help of the previous active contours, as shown in References 53 and 54. Similar results have also been obtained by Herlin et al. with Markov random tields55.

  The advantage is a much faster extraction of the edges directly from the data acquired by an x-ray scanner

Figure 7

 7 Figure 7 Left: MRI cross-section of the body at the level of the heart. Right: segmentation result, obtained by EXCALIBUR, into four regions based on texture analysis (courtesy of Ltvy-V6hel)

  crest lines can be found68s74. The one closer to our definition is probably in the book by Hosaka7j (p. 98).

  Figure15shows the extraction of ridges on the surface of a skull acquired with CT1 in two different positions.The idea of deriving the intensity to extract different characteristic or on isophote surfaces, is also successfully applied by (among others) Harr Romeny et ~l.~', Whitaker72 and Morse et ~1.~~.

Figure 8

 8 Figure 8 Ridges on the brain (courtesy of Thirion and Gourdon)

  Figure IO shows the extraction of ridges obtained from the B-spline approximation.

  and 15 were extracted with the marching line algorithm.

Figure 9

 9 Figure 9 Segmentation of a face with a deformable B-spline within an image of the head acquired with MRI. Original MRI images were provided by Dr Ericke, Siemens (courtesy of A Guiziec)

elFigure 11

 11 Figure 11 Surface ribbon around a space curve: at each curve point, one can compute the curve Frenet frame (1, n, b), where t and n are, respectively, the curve tangent and normal, and the local surface frame (el, ez, N), where e, and Ed are the principal directions of curvatures in the tangent plane of the surface, and N the surface normal (courtesy of A Gutziec)

Figure 14

 14 Figure 14 Recognition: each time a ridge point of the scene is matched with a model point through the previous hash table, the Frenet frames allow the computation of the six parameters of a rigid displacement (R, t). This displacement is stored in a 6-dimension accumulator, which is stored with another Hash table (courtesy of Gutziec)

  Figure 17.

Figure 15 Figure 16

 1516 Figure 15 Top: Ridges (in black) extracted on the surface of a skull scanned in two different positions (left and right). Original 3D images are produced by a GE-CGR CT-Scan (courtesy of Thirion and Gourdon)

Figure 18

 18 Figure 18Left: a surface is extracted from a 3D CT1 image and a curve is extracted from another 2D cross-section.Right: the curve is matched on the surface with the use of differential geometry constraints.Bottom: superposition of the curve in the image plane extracted from the 3D image (courtesy of Gourdon)

Figure 19 ,

 19 Figure 19 A potential-based method allows the registration of an MRI image (top left) with an NM1 image of the head (top right). An interpolated cross-section of the registered 3D MRI is computed, and its edges are shown superimposed on the NM1 images (bottom). Images are courtesy of Jael Travere, from the Cyceron Center in Caen, France (courtesy of Malandain and Rocchisani)

Figure 20

 20 Figure 20 It is possible to use an average representation of the ridge ines to define an 'atlas of ridges', on a skull and on a brain. Then, this atlas can be used to automatically label the most stable lines of an xbitrary skull (left) or brain (right) (courtesy of Subsol and Thirion)

Figure 21 Figure 22

 2122 Figure 21 Inter-patient registration. Left: best rigid; centre: best aftine; right: best locally affme. All these transformations tend to preserve the correspondence between regions of high curvature, and take into account the local modification of the curvature through an affine transformation. By this approach, the result preserves anatomical invariants based on curvature extrema (courtesy of Feldmar)

Figure 23 Figure

 23 Figure 23 From left to right and top to bottom. (a)(b) Segmentation and (c-g) tracking of the mitral valve of the heart in a sequence of US1 using the active physical model of Nastar; (g) point-to-point correspondences obtained. These images were preprocessed with the sonar-space filtering method of Herlin referenced in the text (courtesy of Nastar)

  Figure 25 shows a modal approximation of the 2D motion of the mitral valve shown in Figure 23, while Nastar showed that the motion field of the surface of the 3D ventricle shown in Figure 27 could be approximated with only nine principal modes and five

Figure 25

 25 Figure 25 From left to right and top to bottom. Four high amplitude eigenmodes of the valve (a,b,c,d). Their superposition (e) compared to the expected result (f) (courtesy of Nastarj

.

Figure 27

 27 Figure 27 Spatio-temporal evolution of left ventricle. Left: original ventricle data consist of 18 temporal images of size 98 x 100 x 110 (courtesy of Dr R Robb, Mayo Clinic); centre: synthesized motion with nine principal modes and 18 temporal values per mode; right: synthesized motion with nine principal modes and only five temporal Fourier coefficients (i.e. 45 real numbers) (courtesy of Nastar)

Figure 26

 26 Figure 26 Computation of the motion field computed through a time sequence of 3D CT-scan images of the left ventricle during the systolic stage, using a finite element technique (courtesy of Benayoun). (Original CT images are a courtesy of Dr R Robb, Mayo Clinic)

Figure 31

 31 Figure 31 Augmented reality. Finding the geometric transformation between a 3D object and its projective transformation allows us to superimpose a 3D representation of the brain vessels obtained from MRI (right) with a 2D angiography of the same vessels (left). This is useful to interpret the images during interventional angiography (courtesy of Feldmar. Original images are provided by General Electric Medical Systems Europe, But, France)

  at each ridge point, it is possible to design an extremely efficient matching algorithm.

			This algorithm,	proposed	by
	Gukziec, has two main stages, namely preprocessing
	and recognition,	and combines	geometric	hashing,
	accumulation	and prediction-verifications2	89:
	a<<
	A
	PC<
	A
	7

May 1995

Number 4 May 1995

ACKNOWLEDGEMENTS

Of course, acknowledgements

Medical computer vision, virtual

Figure 28 Top right: a superquadric is initialized to fit grossly the left ventricle observed with nuclear medicine imaging (top left). Free from deformations make the local tit much more accurate at the expense of a few more parameters.

These parameters can then be tracked through time to assess quantitative values useful for the diagnosis (courtesy of Bardinet and Cohen) some complex therapy'7m'9. '39. Among them, craniofacial surgery or laparoscopy, for instance, are desperately seeking for pre-operative simulation627 14'. This can be done with the help of pre-operative 3D image analysis, and the use of advanced interactive 3D graphics involving virtual reality, for instance"' 'L 16.

Currently, in the EPIDAURE project, Delingette and Subsol, with the help of Cotin and Pignon' 14*, are looking closely at the problem of simulating craniofacial surgery, following the ideas of the pioneering work of Terzopoulos and Waters'25. This is illustrated by Figures 29 and 30. Also, Bro-Nielsen is currently

Figure 29 Virtual reality environment for craniofacial surgery simulation (courtesy of Eidelman and Delingette) Figure 30 Top left: an elastic mode1 of the head is designed from the 3D image of the head. Top right: the skull is cut and its shape is modified with a virtual 3D hand. Centre and bottom: the elastic model of head of the patient deforms itself accordingly to the modified shape of the skull, and the expected result (a deformed Herve Delingette!) can be observed (possibly with texture) before any real surgery is done (courtesy of Delingette, Subsol, Cotin and Pignon) trying to model, with adequate tools, the elasticity of human tissues'43.

Another study in the EPIDAURE group is related to augmented reality, acquired by superposing processed preoperative medical images on per-operative medical images (for instance, video images, x-ray images or ultrasound images). Recent results are presented elsewhere14, and also in Report, 1697Report, (1992) )


Chahab Nastar

Nicholas Ayache

Frequency-based Non-rigid Motion Analysis: Application to Four Dimensional Medical Images

Keywords: Medical image analysis, nonrigid motion, deformable models, modal analysis, Fourier analysis, compression, dynamic data, four-dimensional images, cardiac imagery, automatic diagnosis

We present a method for nonrigid motion analysis in time sequences of volume images (4D data). In this method, nonrigid motion of the deforming object contour is dynamically approximated by a physically-based deformable surface. In order to reduce the number of parameters describing the deformation, we make use of a modal analysis which provides a spatial smoothing of the surface. The deformation spectrum, which outlines the main excited modes, can be efficiently used for deformation comparison. Fourier analysis on time signals of the main deformation spectrum components provides a ternporal smoothing of the data. Thus a complex nonrigid deformation is described by only a few parameters: the main excited modes and the main Fourier harmonics. Therefore, 4D data can be analyzed in a very concise manner. The power and robustness of the approach is illustrated by various results on medical data. We believe that our method has important applications in automatic diagnosis of heart diseases and in motion compression.

INTRODUCTION

Motivation and Organization

N this article, we propose a unified approach for non-I rigid motion estimation from time sequences of threedimensional images [5], i.e., 4D data, by taking into account both spatial and temporal frequencies of a deformable geometric model. Our method has important applications in automatic diagnosis of heart diseases and in 4D data compression, as we shall see in the experimental section. Our method involves three steps:

I) Recover the deformation field between each pair of 2) Express the modal coefficients (or amplitudes) of the 3) Express Fourier coefficients of the time-varying mo-

The information contained in the data can then be compressed, by discarding spatial modes and/or temporal Fourier harmonics.

The basic justification of the approach relies on the following observations: most smoothly deforming structures mainly have low-frequency excited modes; this justifies Step 2. Furthermore, for periodic motions like heart motion, modal amplitudes as a function of time are periodic and sine-like; this justifies Step 3. successive 3D images, deformation at each time t, dal amplitudes.

For Steps 1 and 2, it is possible to constrain the modes beforehand, by letting a deformable model evolve along loworder modes. The advantage is, of course, a reduced numerical complexity and a low-order smoothing of the deformation allowing its robust recovery. Note that Step 1 can be performed not only by our deformable model, but also by optical flow or any other technique providing the motion field [lo].

The following example can serve as a motivation to our work. Let us consider (Fig. I) the canine 4D heart data provided by the dynamic spatial reconstructor (DSR, a high speed X-ray CT scanner). We used this data as an input to our method. It consists of 18 volume (or 3D) images during a single cardiac cycle. Each volume image has a spatial resolution of 98 x 100 x 110. Dye was injected into the left ventricle, which shows up as a light gray color. Note the difficulty of interpretation of such complex and huge data by physicians. The challenge was to analyze such a tremendous amount of information [START_REF] Danielsson | Euclidean Distance Mapping[END_REF]404,000 bytes, with each voxel being coded on one byte) and supply physicians with a few quantitative parameters describing the motion. We shall see in the experimental section of this paper how we tried to meet these objectives.

The article is organized as follows:

In Section 2, we introduce our deformable model for nonrigid motion estimation, develop its governing equations and outline its properties.

In Section3, we present the modal analysis of the model, allowing a closed-form recovery of the motion by a few parameters: the modal amplitudes. We show how these amplitudes may be used to characterize deformations. Section 3.2, which describes the analytic modes, may be required only for in depth reading.

(step 3).

Related Work

A lot of research work is related to our work. First of all, we want to mention the pioneering work of Terzopoulos et al. on snakes and deformable physics-based models [241, 1261, [42], [START_REF] Terzopoulos | Constraints on Deformable Models: Recovering 3D Shape and Nonrigid Motion[END_REF]. This work was mainly dedicated to tracking and animation problems rather than analysis.

A second category of research work is related to the analysis of deformations. Our work was initially inspired by the pioneering work of Pentland et al. on the use of modal deformations to describe deformations [START_REF] Pentland | Recovery of Non-Rigid Motion and Structure[END_REF], [371. The main differences with our approach is that our modes have a simpler physical meaning (they can be interpreted as the harmonics of the free vibrations basis of the deformable object) while Pentland's modes are more dedicated to animation [38]. Also, we have derived a method to compute our modes beforehand, as a function of the topology of the deformable object only, allowing a much faster computation (Pentland's modes require the knowledge of the exact shape of the original object and, therefore, have to be computed online).

Our work can be seen as providing a coherent unifying framework between these first two categories of publications, allowing to both track and analyze a deformable 3D motion.

In addition to these references, we have to mention the work of Cootes et al. [15], [161, [17], done independently of, and concurrently with, ours which provides an interesting alternative to the solution of both the tracking and the analysis of the deformation. We also have to mention the work of Bardinet [7] and Metaxas [27], [START_REF] Park | Volumetric Deformable Models with Parametric Functions: A New Approach to the 3D Motion Analysis of the LV from MRI-SPAMM[END_REF], among others, who use deformable parametric models like superquadrics to track and analyze deformable objects. Their work is limited to the class of objects which can be easily represented by such parametric models, which include the left cardiac ventricle. A number of methods to track 3D data were proposed in the past few years. Duncan et al. 121,[START_REF] Shi | A Unified Framework to Assess Myocardial Function from 4D Images[END_REF], 1401, and Benayoun et al. 191,[lo], used differential constraints; Goldgof [23] and Chen 1131 used a coarse-to-fine approach; Bookstein [12] and Szeliski and Lavallee [START_REF] Szeliski | Matching 3D Anatomical Surfaces with Non-Rigid Volumetric Deformations[END_REF] proposed energy based methods; Creswell [START_REF] Creswell | Mathematical Modeling of the Heart Using Magnetic Resonance Imaging[END_REF] and Amini [l] suggested data-dedicated methods. All these approaches can be seen as complementary to ours, in the sense that the tracking result they provide could be completed by our proposed analysis. It would however be interesting to compare all these approaches on a specific and common database. This could be the topic of a future work. In Section 4, we propose a Fourier analysis for the timevarying main modal amplitudes; this is particularly well suited for cyclic motions (e.g., heart motion). An example involving 4D medical data is presented.

We conclude and propose future work in Section 5.

A DEFORMABLE MODEL FOR MOTION ESTIMATION

In this section, we introduce our physically based deformable model which we use for tracking nonrigid motion of dynamic structures in time sequences of 2D or 3D medical images. We consider both the surface and volume properties of the objects at hand. We restrict ourselves to elastic deformations, i.e., we assume that the object recovers its reference configuration as soon as all applied forces causing deformation are removed. In general, we seek a trade off between precise modeling and computational efficiency. Therefore, simplifying assumptions will be introduced in the modeling.

Mass-Spring Meshes

-4

where K is the stiffness constant, U,(t) = MlMl(t) is the vector separation of nodes M , and M , at time t, and C, Modeling an elastic boundary can be achieved by a mesh of N virtual masses on the contour. Each mass is attached to its neighbors by perfect identical springs of stiffness K and natural length lo [231, [461. These springs achieve a polygonal approximation and I is the set of nodes connected to node Mi.

Fluid damping force, proportional to node velocity: model the surface properties of the object.

Generalizing the model to 3D contours (surfaces) is straightforward: We can either model "quadrilateral" or "diagonal" meshe? (see Fig. 2).

Fig. 2. "Diagonal" mesh

If necessary, we can improve the modeling by attaching extra springs between non-neighbor nodes in order to model volume properties inside the object. These springs constrain the general shape of the object within its deformation, making a template shape out of it [ll], [481.

The boundary modeled as above will also be called a stuucture. Such a structure can be easily deformed to match the contour of an object of interest, thus performing a segmentation step. If we take a series of images displaying the deformation of the object through time, the structure can achieve simultaneously both segmentation and tracking of the object surface through time.

Governing Equations

The system under study is made of the N virtual masses located at time t at points (Ml(t), M2(t), . . ., MN(t)).

The fundamental equation of dynamics states that the vector addition of all applied forces on M, is equal to its mass m, multiplied by its acceleration. Let A be the origin of the reference system. Then:

Let us now describe the applied forces on each node Mi:

Elastic force between M i and its connected nodes.

(2)

f,(M,,t) = -C -AM! d -3 = -c,M, ' dt (3)
where c, is the damping constant.

Image force flm (Mf, t ) . This is the main external force that has the system attracted by image features. This force is defined in Section 2.3. Suppose the natural length 1, of the springs is fixed. Since we wish to give the system an initial equilibrium configuration, we need to apply on each node a force balancing the action of the elastic force: f,,(MJ = -fe(Mz, to)

(4)

This force is similar to the force that our fingers apply to an elastic rubber to keep it in a specific shape. We assume that this force is constant over time.

Finally, (1) yields the governing equation:

fc(Mr, t ) + fd(Ml, t ) + & ( M I , t) + f,(M,) = m,M, (5)
The governing equation, expressed for all N nodes, leads to a nonlinear system of coupled differential equations (for each node, the x, y, and z displacements are coupled, and the displacement of a node depends on its neighbors displacement, as it appears clearly in (2)).

One possible approach is the resolution of these complex equations by an iterative procedure [43]. In this paper, we propose to set 1, = 0. This assumption does not restrict the arbitrary initial configuration of the structure because of the equilibrium force feq. Indeed, this force keeps the structure inflated so that it does not shrink to a point. Thus, the natural state of the system is its initial configuration.

The advantage of this assumption is that our model can be considered within the framework of linear elasticity. As a consequence, we end up with a set of linear differential equations with node displacements decoupled in each coordinate, regardless of the magnitude of the displacements. Moreover, these linear equations are a prerequisite to further quantitative analysis of the motion (see Section 3).

On the other hand, our approximation is valid only if the spring orientations undergo small angular variations (typically less than 15 degrees), so that our assumption of constant equilibrium force f,,, holds.3

Finally, in 3D, the deformation of the system is governed by the 3N-dimensional differential matrix equation:

MU + CU + KU = F(t) (6)
where U is a vector storing nodal displacements M, C, and K are, respectively, the mass, damping, and stiffness matrices of the system, and F is the image force which has the object attracted by image edges. Equation ( 6) is the finite element formulation of the deformation process. Note that the equilibrium forces do not explicitely appear in the governing equation.

.

The shear-resisting cross springs may be useful to avoid selfintersection problems during the deformation of the surface.

3. Similar limitations can be found in the model described in [361, 1371

Image Force

In the original formulation of the "snake," the authors propose a convolution of the image with a smoothing filter that causes artificial blurring, so that the active contour can be attracted by the edges from a distance [24]. Unlike these methods, we introduce a force at each node [START_REF] Ye | The Signed Euclidean Distance Transform and Its Applications[END_REF]. At node M, this force is set to:

, (7) fi,(M, t ) = (X M(t)P(t)
where a is a constant scalar. This force can be seen as a virtual spring of natural length zero and of stiffness a joining M to P. Hence, both internal and external forces are elastic, making the modeling coherent.

The advantage is that we speed up the convergence of the model toward image edges; one can consider edge extraction and distance computation as a data-to-force transformation which is performed as a preprocessing.

Integration Scheme

In 3D, the 3N-order matrix equations decouple into three directional matrix equations of order N:

MU, + CU, + KUx = F,(t) MUz + CUz + KUz = F,(t) MU, + CU, + KU, = F,(t) (8)
larger time step can be used when high-frequency components of the deformable model are discarded (see Section 3).

Results on Medical Data

We have tested our method on a set of ultrasound images of the left ventricle of a human heart. The tracking of the mitral valve is indeed a problem of major interest in medical imaging since heart-attacks can generally be predicted from abnormal motion of the valve. Each image has a resolution of 256 x 256 pixels. First, a polar edge extraction is performed on the images [61. Then, for each image, the distance field is computed on every pixel. Thus we can segment the initial valve. We then add volume springs for making our model a template; this will make the tracking more robust. Finally we track the valve through time and display the estimated displacement field on the valve surface (Fig. 3). The program runs in real-time on 2D data.

where M, C, and K are from now on N-order matrices, and Up and Fp (p = x, y, z) are, respectively, the N-order displacement field and image force in the ,u direction. From now on, we will omit the indexes, and all matrix equations will be of order N. It is assumed that three equations corresponding the three space directions have to be solved.

We integrate the governing equations with an explicit Euler scheme:

U(t) = M-'(F(t) -CU(t) -KU(t)) U(t + At) = U(t) + AtU(t + Af) U(t + At) = U(t) + AfU(t) (9)

I

where At is the time step of the simulation. The initial values of displacement and velocity are generally set to zero. Note that the mass matrix is diagonal and therefore its inversion is trivial and only has to be processed once. The damping constant is chosen so that the system is slightly overdamped [44]. Convergence is achieved once the image force balances the internal elastic forces.

In any numerical integration scheme, numerical stability requires the time step At to be inferior to a critical value Atcr which is defined by the mass and stiffness properties of the system:

Y

where TN = * is the smallest period of the finite element mesh, and yis a constant depending on the scheme. The Euler method is conditionally stable, i.e., needs small time steps ( y 2 2).

order to make a more accurate computation, we decrease the time step and set y 2 10. Note that a ai The method has similarly been tested on a 3D Magnetic Resonance image of the human head, with a resolution of 158 x 158 x 158 voxels. We make use of the 3D generalization of Canny's edge detector [28], 1291, then compute the 3D distance maps and finally run the program. The convergence of an initial sphere toward the human head is shown in Fig. 4. It runs at interactive rates on a DEC-ALPHA workstation; convergence is achieved in eight seconds. Note that this segmentation example does not really make use of the power of our dynamic approach; it is provided here to demonstrate model properties.

MODAL ANALYSIS AND APPLICATIONS

In this section we focus on a quantitative analysis of the model deformation, using a frequency-based technique called modal analysis.

Modal analysis is a standard engineering technique allowing more effective computations and a closed-form so- lution of the deformation process [8]. It was first introduced in computer vision by Pentland's team [361, 1371. Let us now explain this technique.

General Approach

Instead of solving directly the equilibrium equation ( 6), one can transform it by a change of basis:

U = l % (11)
where P is the square nonsingular transformation matrix of order N to be determined, and lb is referred to as the generalized displacements vector. One effective way of choosing P is setting it to @, a matrix whose entries are the eigenvectors of the generalized eigenproblenc KQ= w2MQ

(12

) N u(t) = @tj = C.",(t)$b, ( 13
) I=1
Equation ( 13) is referred to as the modal superposition equation. is the ith mode, iiiits amplitude, and q its frequency. Fig. 5 displays a sample of six frequency-increasing modes of a cylinder, all of them having the same amplitude. Note that the modes of a generalized eigenproblem involving real symmetric matrices can be chosen to be orthonormal vectors [8]. The new modal basis simultaneously diagonalizes M and K, and provided that matrix = QTC@ is diagonal as well,4 the governing matrix-form equations decouple into N scalar equations:

ii,(t) + q g t ) + w: u",(t) = I ( t )

(14)
The amplitudes (Zl(t))lzl, ,N, are obtained by solving these equations at time t, and the displacement of the structure nodes is obtained by the modal superposition equation.

In practice, we wish to approximate nodal displacements U(t) by u ( t ) , the truncated addition of the p low-frequency modes, where p < N.

(15)

Vectors (@Jixl,. ..,i) form the ueduced modal basis of the system. This is the major advantage of modal analysis: It allows a closed-form solution by selecting a few number of low-frequency modes [8], 1311, 1321. Therefore a compact description of the motion is provided by spatial smoothing. Fig. 6 illustrates modal superposition on the real example of the mitral valve.

4. This condition, called the Rayleigh condition, is satisfied as soon as the damping matrix C is a linear combination of the mass and stiffness matrices. Rayleigh damping is generally assumed for any standard engineering problem 181.

initially a sphere, one can observe that the error is more important in the regions that need a locally important deformation of the sphere (nose, eyes, ears), while more spherical regions (top of the head) are better approximated, as expected. Note however that the level of accuracy can be easily controlled, at the expense of additional modes.

Analytic Modal Analysis

Fig, 6. Nine low-order amplitude-moderated eigenmodes describing the motion of the valve; their superposition (bottom-right) is a good approximation of valve motion.

Note that superimposing p low-frequency equations of the type ( 14) allows the time step to be larger; we choose: ?r 2 z whereT = -B T, . As a consequence, the convergence of the model is numerically much faster with modal analysis than with direct integration of the governing equations.

Due to orthonormality of the modes, the approximation error is:

(1 7 )
This means that the approximation error is a rapidly decreasing function of the truncation frequency p. (Exact computation needs 37,683 modes).

The distribution of the approximation error along the surface for increasing values of p is shown in Fig. 7. Darker colors outline regions with bigger error. The model being Even as a precalculation, solving the generalized eigenproblem is very costly as soon as we consider 3D boundaries (surfaces). For instance, if we consider a mesh of 100 x 100 nodes, a generalized eigenproblem where the size of the matrices is 10,000 x 10,000 has to be solved. It is clear that the analytic expvession of the modes would noticeably reduce the computations [30]. This leads us to consider the theory of solid state physics, where similar problems are encountered at a microscopic level (ionic vibrations of a crystal lattice). If we parameterize our deformable curves by arc length, and similarly our deformable surfaces by natural coordinates, we get periodic boundary conditions which depend on the surface topology. This allows the analogy between our deformable model and a crystal lattice.

Free Vibrations of a Chain

The classical theory of vibration of a crystal lattice is based on the harmonic appuoximation, a theory which assumes that the first non-vanishing correction to the equilibrium potential energy is quadratic [3]:

(18) 1 2 _ . Vharm = -x u (l')Duv(Z -i?')pv(Zr) R,R' U,V
where ~~( 2 ) is the displacement in the , U direction of the ion whose mean position is R , and D is the Hessian matrix of the interaction energy.

Consider now a set of ions distributed along a closed chain at points separated by a distance a, so that the lattice vectors are R = nu for n E (1, . . ., NI. If only neighboring ions interact, we may take the harmonic potential energy to have the form: n=1 where K = v"(a) is the stiffness constant of the system, v(x) being the interaction energy of two ions at distance x along the chain. The free vibrations of the lattice are governed by: (20)

dV knvm Mii(na) = -___ ~ -K(u((n + 1)a) + U ( ( . -1)a) -2u(na))
These are precisely the equations that would be satisfied if each ion were connected to its neighbors by perfect massless springs of stiffness K (and equilibrium length a, although the equations are in fact independent of the equilibrium length of the spring).

We seek solutions to (20) of the form Ae'ikna-wi). This yields the dispersion equation which gives the relationship between spatial ( k ) and temporal ( U ) frequencies: The periodicity of the closed chain is expressed by u[(n + N)al = u(na). We now obtain the set of independent solutions:

k(p)a = __ 'r p E B ( N ) (22)
where B(N) is the first Brillouin zone. B(N) is equal to {-+ t 1 , ...,+} for N even, and

The general solution of the free vibrations of the closed {p , ..., y } for N odd. chain is the linear combination of the former solutions:

(23) (24)
The case of the open chain is very similar. We sum up the expression of the free vibrations for both types of chains in , _ Table 1.

TABLE 1 DISPLACEMENT OF THE NODES OF A FREE VIBRATING CHAIN

Nonlinear Waves in Discrete Media

When k is small compared with n/a (i.e., when the wavelength is large compared to the interparticle spacing), w is linear in k (from (21)):

w = & (25)
This is the type of behavior we are accustomed to in the case of light waves and ordinary sound waves. If w is linear in k, then the group velocity is the same as the phase velocity (equal to c = a g ) , and both are independent of frequency. Note that if we approximate finite differences by derivatives:

u((n + 1)a) -u(na) = adu/dx(na) u((n -1)a) -u(na) = adu/dx((n -1)a) u ( ( n + 1)a) -u(na) + u((n -1)a) -tl(na) = a2~2ti/ax2
in (20), we end up with a wave equation of velocity c: [START_REF] Mcinerney | A Dynamic Finite Element Surface Model for Segmentation and Tracking in Multidimensional Images with Application to Cardiac 4D Image Analysis[END_REF] -

with c = a - F d 2 U 1 d2u ax2 c2 at2 --_
One of the characteristic features of waves in discrete media, however, is that the nonlinearity ceases to hold at wavelengths short enough to be comparable with the interparticle spacing. In the present case w falls below ck as k increases, and the group velocity drops to zero when I k I reaches d a .

Analytic Modes for Curves

In the more general case of damped and forced vibrations (which is the case of our governing equations), time dependency is not harmonic, andl has to be computed separately. For a closed chain, (23) becomes:

u ( m , f) = c i i , ( t ) P P 4 N ) (28)
Comparing ( 28) with the modal superposition (13) yields the analytic expression of the modes for a closed chain. Note that since time and space dependency are separate in (281, modal analysis is the decomposition of the displacement in a basis of standing waues which are the vibration modes of the system.

Analytic Modes for Surfaces

The generalization to surface meshes is done by mixing all possible pairs of boundary conditions; it yields three different surface topologies: openandopen

+ torus topology, closed and closed plane topology, closed and open

+ cylinder topology

The analytic expressions for these topologies are summed up in Table 3, where the set of variation of the mode parameters p and p' is %(N) or (0, ..., N -11, depending on the boundary conditions.

Note that we implicitly develop the surface expressions for quadrilateral meshes, and that the modes have to be normalized to unity.

Finally, analytic modal analysis has both theoretical and practical implications. Theoretically, it shows that modal analysis is a specific form of Fourier decomposition of the deformation in a basis of standing waves (28). From a more practical point of view, analytic expressions of the modes are an efficient tool for veal time eigenvector extraction as soon as the surface topology and the mass and stiffness properties of the model are defined.

Analytic Modes for Volumes

Suppose we wish to model a deformable volume. This may be typically the case for nonrigid motion recovery in a whole 3D image [lo]. Therefore we model a volume mesh of size N x N' x N". The volume has "open" boundary conditions in all three directions. Thus, it is easy to generalize the results for the deformable plane, defining this time the modes with three pa-plane torus cylinder deformation spectrum as the set of total amplitudes of the first p deformable and low-frequency modes. Once the spectra are computed, we can define a distance measure between the spectra. We choose the Euclidean distance d, such that lower

The eigenvalues are:

2N

(29) amplitudes are given less importance than higher ones:

while the eigenvectors have the following expression:

4(Pr P'r P") = [START_REF] Nastar | Classification of Nonrigid Motion in 3D Images Using Physics-Based Vibration Analysis[END_REF] where D, and D2 are the labels of the deformations. Dis- [START_REF] Nastar | Vibration Modes for Nonrigid Motion Analysis in 3D Images[END_REF] tance d gives a relative value of how different the two de-PZ(2fl-1)

(
P'Z(2n' -l) P"Z(2n" ~ 1) cos 2N" , . . . ]
2" formations are. [...,cos 2N cos

with n E (1, ..., NI, n' E (1, ..., MI, and n" E (1, ..., N"].

I Similarity

Fig. 8 shows two similar deformations and their spectrum; in this example, R is a 90-degree rotation matrix and T is 3.3 The Deformation Spectrum and Its Applications DEFINITION. The deformation spectrum of a motion is the graph representing the value of the modal amplitudes as a function of mode rank: ct(t) = f(i) [START_REF] Nastar | Classification of Nonrigid Motion in 3D Images Using Physics-Based Vibration Analysis[END_REF]. The deformation spec-arbitrary.

model for segmenting the mitral valve contour in Fig. 3.

These are the deformations of the initial deformable trum is initially drawn for a deformation occurring between two image frames: It describes which modes are excited, and how, in order to deform one object into another.

It also gives an indication of the strain energy, as we have:

1 1 N (31) Estrain = -U'KU = -2 c4.",' !=1 2
Note that rigid motion has zero strain energy. Let us now define the term similar deformations. Two deformations are similar when the corresponding displacement fields U, a rotation matrix R and a translation vector T such that:

(1 U1 -(RU, + 73 ( 1 < E and U, are similar within a rigid transform, i.e., if we can find

Reference lMormsllon for a small value of E.

Analytic Therefore, it is natural to state that, provided that the dimensionality p of the reduced modal basis is suitably chosen, two similar deformations have similar deformation spectra. However, we cannot compare the deformations of two objects placed in arbitrary configurations in the same global reference frame. Thus the modal computations have to be developed in the object reference frame, defined by its center and axes of inertia.

In the following, we get rid of the rigid modes in order to study exclusively the deformations. Thus we will refer to the

Robustness

In order to test the robustness of the spectra in the presence of noise, we add a Gaussian noise of standard deviation o to the initial 3D diastole, then deform the diastole into the systole, draw the spectrum and compare it to the spectrum of the true deformation.

The results are very promising (Fig. 9): For o = 0.1 the difference between the original spectrum and the corrupted one is almost invisible (d = 0.01); for g = 1 very slight differences appear (d = 0.10); for G = 5 the main excited modes are still the same (d = 0.44).

Note the very chaotic visual representations of the corrupted shapes, while the modal spectra still succeed in extracting the significant deformation information.

Classification

In order to test the matching of a specific deformation to a group of predefined deformations, we can consider several admissible deformation? X, of a reference shape like the valve (Fig. 10). Under the assumption of Gaussian distribution, we can then classify a test deformation Y as belonging or not to the set D = {. . ., X,, . . . J by using the Makalanobis distance [41,[211:

dL(Y,'Z)) = (Y -X)'W-l(Y -x) (34)
where W is the covariance matrix of the admissible deformations D, and x the mean admissible deformation.

By comparing with a x table of q degrees of freedom, where 9 is the rank of the covariance matrix,' we can determine a confidence measure for acceptation

( d i ( Y , D ) < E )
or rejection of the test deformation as being part of the predefined admissible deformations. Fig. 11 shows the classification of four deformations. The confidence that we have in our classification (derived from the x ' table) is also indicated. Note that visually, we would classify the four deformations the way the system has done it (that is, reject the first three deformations, and accept the fourth). Note finally that though this set of deformations has been artificially generated, the application of the method to a clinicallysignificant case is straightforward.

2

TIME EVOLUTION OF THlE MAIN MODAL AMPLITUDES

A single deformation spectrum gives a static information about the spatial frequencies of the motion, whereas the time parameter t is not really taken into account.

However, temporal evolution is really what we are interested in. A dynamic deformation process can be much better described and interpreted if we have a time sequence of images. Therefore, we can draw T consecutive deformation spectra for T + 1 frames of images showing the temporal evolution of the deformation process. For a chosen mode i, we are interested in the time signals: Gl(t) for different values of i.

5. By deformation we mean the vector of total modal amplitudes, in other 6. If 9 < p then W -1 is the pscudo-inverse of W.

terms the deformation spectrum. Let us track the moving ventricle (cf. Fig. 1) through all image frames. We wish to find a minimum number of parameters describing this sequence. Can we compress the 4D data information for its further analysis, storage, or transmission?

Let us denote by (Po, . . ., Pr) the positions of this surface during the T + 1 frames of the cardiac cycle. Therefore, 3 N x (T + 1) parameters describe the 4D data (i.e., 345,600 parameters with N = 6,400 and T = 17).

Let us now track the ventricle in the reduced modal space. We obtain T truncated spectra that store the lowfrequency modal amplitudes through the cycle. 4D data is now stored in 3N + T x p parameters (i.e., 20,985 parameters with p = 105). The numerical value of p is chosen so that each truncated spectrum has 90% of the energy of the corresponding non-truncated spectrum.

We wish to discard as many modes as possible. Therefore we define a criterion for keeping only the most excited modes among the p low-frequency modes. We compute the energy of an amplitude through the cycle: T L2(i) = I, ;r,(t)*dt [START_REF] Park | Volumetric Deformable Models with Parametric Functions: A New Approach to the 3D Motion Analysis of the LV from MRI-SPAMM[END_REF] The larger L2, the more important the contribution of the corresponding amplitude. In our experiment with this data, the first nine values of L2 were much larger than the other ones. q = 9 was then an obvious threshold for us to choose.

Let us draw the temporal evolution of the q selected modal amplitudes. They are represented by groups of three in Fig. 12. These curves have, globally, a single period sine shape, which is an expectable result for low-frequency modes during a cardiac cycle. Their shape is quite like the ventricle volume curve as a function of time during a cardiac cycle. The shape of these curves encourage us to perform a fast Fourier Transform: Most probably only a few Fourier harmonics will describe the time evolution of the curves.

, m a -V In order to illustrate this point, let us choose one among these q = 9 spectral components (or modal amplitudes). As for motion compression, instead of keeping a whole 3D image at every frame (that makes t, t, t, =1,078,000 parameters per frame) we now keep only H x q/(T + l) = 2.5 parameters per frame. This compression of 4D data is impressive: It means that, provided that we keep the first shape, we are able to synthesize the motion with 2.5 parameters per frame. This 4D data is then described by 3N + H x q = 19,245 parameters instead of (T + 1) x t, t t = y 2 3

19,404,000, which makes a compression of more than 1 0 , for a cardiac cycle whose duration is approximately one second. If the analysis was extended over a much longer period (e.g., 10 minutes), the Compression rate would continuously increase toward the asymptotic value of 1,078,000/2.5, i.e., approximately 4.10 per pame. If we compare the total compression with the transmission of all the 3D images over the 10 minutes, we would get an asymptotic compression factor superior to 10 .

On the top of the compression ability, the synthetic information that our spatiotemporal processing provides allows comparison of nonrigid motion by comparing very few parameters. The method is indeed useful not only for compression, but also for analysis of dynamic motion for diagnosis purposes.

The evolution of the left ventricle model over a cardiac cycle is displayed in Fig. 15. In each row, we observe: Left, the evolution of the model in the real space (6); center, the approximation with 9 = 9 modes (spatial smoothing); right, spatiotemporal processing with nine modes and five Fourier harmonics.

Finally, Fig. 16 displays the mean Euclidean error between the real-space mesh (complete computation) and respectively spatial smoothing (modal recovery) and spatiotemporal smoothing (modal and Fourier recovery). Note that this error (in voxels) is extremely small. The position of the two peaks indicate that the error is maximum in the middle of diastolic and systolic shapes: This is because the change in the overall shape is maximum in those frames. An interesting topic would be the study of the error range, and its relationship with the number of modes and Fourier harmonics.

CONCLUSION

In this article, we presented a frequency-based analysis of nonrigid motion by coupling modal analysis (outlining the main spatial frequencies of the motion), and Fourier analysis (outlining the main temporal frequencies, mainly for cyclic motions).

Nonrigid motion is first estimated by a physically based deformable model.

The modal decomposition of the deformation, made real-time and precise by using the analytic expressions of the modes, leads to the definition of the deformation spectrum, which is a compact description of the deformation allowing its straightforward comparison and classification (identifying pathologic and normal deformations).

We then introduced a temporal analysis of nonrigid motion from 4D data by recovering the temporal evolution of the main modal amplitudes by a Fourier analysis.

The method provides a very nice description of 4D data by very few parameters (a few modal amplitudes and a few Fourier harmonics). It has important applications in medical analyis of nonrigid motion (mainly for automatic diagnosis purposes), and also provides a tremendous compression of multidimensional nonrigid motion, mainly for storage or transmission purposes.

Our future work will focus on the clinical use and validation of this framework, and its transfer toward concrete medical applications, in particular in the cardiology field.

Fig. 1 .

 1 Fig. 1. 4D data displaying a 3D dog ventricle in motion. Each 2D image is an axial slice in the x y plane. The z coordinate is along the vertical axis, while the tcoordinate is along the horizontal axis.

Fig. 3 .

 3 Fig. 3. Initial segmentation and tracking of the mitral valve

Fig. 5 .

 5 Fig. 5. A sample of frequency-increasing modes of a cylinder (constant amplitude).

Fig. 7 .

 7 Fig. 7. Approximation error distribution for the low-order deformation of a sphere into the human head. Darker regions indicate more important errors, while lighter regions represent best-recovered ones. a. Computation using 60 modes. b. 243 modes. c. 468 modes. d. 2790 modes. (Exact computation needs 37,683 modes).

Fig. 8 .

 8 Fig. 8. Similar deformations D, and D2 and their spectrum.

Fig. 9 .

 9 Fig. 9. Corrupted initial diastole deformed to the systole (not shown here), and its spectrum, for different values of noise D(G= 0.1, o= 1, o= 5).

Fig. 10 .Fig. 11 .

 1011 Fig. 10. Ten admissible deformations of the valve.

Fig. 12 (Fig. 13 .Fig. 12 .

 121312 Fig. 12 (bottom right) shows the time evolution of this particular amplitude. Fig. 13 displays the real (top left) and imaginary (top right) parts of the corresponding Fourier Transform; we observe that indeed only low-frequency Fourier harmonics are excited. Therefore we keep the harmonics of rank 0, 1, and 2 of the Fourier spectra (thus, with symmetry considerations, a total of H = 5 harmonics, see Fig. 13 bottom left and bottom right); we then reconstruct the time signal in Fig. 14. nul p r t 1m.Zlmry par,

Fig. 15 .

 15 Fig. 15. Space-time evolution of the left ventricle model. In each row: Left: Finite element formulation. Center: Modal recovery. Right: Modal plus Fourier recovery.

X

  Fig. 16. Evolution of the approximation error (in voxels) during a cardiac cycle, for modal recovery (solid line) and modal plus Fourier recovery (dotted line).

  Table 2 sums up the frequencies (eigenvalues) and modes (eigenvectors) for closed and open curves.

		TABLE: 2
	ANALYTIC EIGENVALUES AND EIGENVECTORS
	FOR CLOSED AND OPEN CURVES
	of x M/4K
	closed	[ ..., c o s y , . . . ' 7
	open	[...,cos	px(2n -1) 2N

I ...p

  and p" varying, respectively, in (0 . . . N -11, (0 ... N'-l}, and (0 ...

		TABLE 3
		ANALYTIC EIGENVALUES AND EIGENVECTORS
		FOR DIFFERENT SURFACE TOPOLOGIES
		U' x M / M	6P)
	2 PZ sin -+sin 2N	2 P' " 2N	pn(2n -1) [...,cos 2N	cos	p'z(2n' -1) 2 "	T
	2 P" sin -+sin N	2 P'Z N
	2 P E sin -+sin 2N	2 P' " N	p 4 2 n ~ 1) [ ..., cos 2N	2p' "
	rameters p, p',

N" -1).

TABLE 3

 3

	.L3ms-
	' I ' ~ ~
	, : -	_ya/	,lla/	11/10	YMX,

62-8828/96$05.00 01996 IEEE

The authovs aye with INRIA, B.P. 105,78153 Le Chesnay Cider, France E-mail: {chnhab.nasfar; nicholas.ayache)@inria.fr. Manuscript received Mar. 14,1995; revised July 25,1995. Recummended for acceptance by A. Singh. For infurriation on obtaining reprints of this article, please send e-mail to: transpami@computer.org, and reference IEEECS Log Number P96085.

ACKNOWLEDGMENTS

We wish to thank Dr. Richard Robb at Biomedical Imaging Resource, Mayo Foundation/Clinic, for providing the DSR data. This work was supported in part by a grant from Digital Equipment Corporation.