H Delingette

X Pennec

L Soler

J Marescaux

N Ayache

Computational Models for Image Guided, Robot-Assisted and Simulated Medical Interventions

Keywords: Computational Models, Medical Robotics, Medical Image Analysis, Anatomical Models, Biomechanics, Physiological Models, Brain-shift, Augmented Reality, surgery simulation Pathology Muscles Activation Respiration Physiological Modeling Statistical Analysis Physical Modeling kinematics Temperature Deformation Electro-magnetism Liquid/Solid

Medical Image Analysis plays a crucial role in the diagnosis, planning, control and follow-up of therapy. To be combined efficiently with medical robotics, Medical Image Analysis can be supported by the development of specific computational models of the human body operating at various levels. We describe in this article a hierarchy of these computational models, including the geometrical, physical and physiological levels, and illustrate their potential use in a number of advanced medical applications including image guided, robot-assisted and simulated medical interventions. We conclude this article with scientific perspectives.

I. INTRODUCTION

Medical Imaging can be used for the planning and control of the motion of medical robots. Examples include the use of pre-operative images and geometric reasoning for path planning, the combined use of pre-operative and intra-operative images to control a surgical procedure with augmented reality visualization, or the simulation of surgical interventions on virtual organs built from pre-operative images and atlases with realistic visual and force feedback.

Such procedures require the use of advanced medical image analysis methods and the development of a hierarchy of socalled computational models of the human body [START_REF] Ayache | Computational Models for the Human Body, ser. Handbook of Numerical Analysis[END_REF]. These computational models aim at reproducing the geometrical, physical and physiological properties of human organs and systems at various scales (see Figure 1). They can be used in conjunction with medical images and robotics to actually enhance the possibilities of image analysis and robot control.

The purpose of this article to describe first a hierarchy of levels of these computational models, and then to illustrate their potential use for medical robotics applications.

The first level is mainly geometrical and addresses the construction of digital static descriptions of the anatomy, often based on medical imagery. The techniques for segmenting and reconstructing anatomical and pathological structures from

Statistical Analysis

Fig. 1. Hierarchy of computational models of the human body (from [START_REF] Delingette | Soft tissue modeling for surgery simulation[END_REF]) medical images have been developed for the past 15 years and have brought many advances in several medical fields including computer-aided diagnosis, therapy planning, image-guided interventions, drug delivery, etc. An distinctive achievement of computational anatomy has been carried out by the "Visible Human Project" [START_REF] Ackerman | The visible human project[END_REF] which provided the first digital multimodal anatomical representation of the full human body.

A second level of modeling describes the physical properties of the human body, involving for instance the biomechanical behavior of various tissues, organs, vessels, muscles or bone structures [START_REF] Haug | Human models for crash and impact simulation[END_REF].

A third level of modeling describes the functions of the major biological systems [START_REF] Hunter | Integration from proteins to organs: the Physiome project[END_REF], [START_REF] Noble | Modeling the Heart, from genes to cells to the whole organ[END_REF] (e.g. cardiovascular [START_REF] Belik | Computational methods for cardiac electrophysiology[END_REF], [START_REF] Sermesant | An electromechanical model of the heart for image analysis and simulation[END_REF], respiratory [START_REF] Kaye | A 3D virtual environment for modeling mechanical cardiopulmonary interactions[END_REF], digestive, hormonal, muscular, central or peripheral nervous system, etc.) or some pathological metabolism (e.g. evolution of inflammatory or cancerous lesions [START_REF] Swanson | Virtual and real brain tumors: using mathematical modeling to quantify glioma growth and invasion[END_REF], formation of vessel stenoses [START_REF] Quarteroni | Mathematical modeling and numerical simulation of the cardiovascular system[END_REF], [START_REF] Boissonnat | From arteriographies to computational flow in saccular aneurisms: the INRIA experience[END_REF], etc.). Such physiological models often include reactive mechanisms while physical models only provide a passive description of tissues and structures. A fourth level not depicted in Figure 1 would be cognitive, modeling the higher functions of the human brain.

Each model is specified by a number of parameters and a related task consists in studying the variability of those parameters across a given population. Thus, statistical modeling of those computational models can be seen as an orthogonal modeling activity that aims for instance at finding the local or global similarities and dissimilarities of a structure or a function between two populations [START_REF] Thompson | Guest Editorial[END_REF], [START_REF] Miller | Computational anatomy: shape, growth, and atrophy comparison via diffeomorphisms[END_REF], [START_REF] Csernansky | Computational anatomy and neuropsychiatric disease: probabilistic assessment of variation and statistical inference of group difference, hemispheric asymmetry, and time-dependent change[END_REF]. Statistical findings may also be used to calibrate, refine or constrain [START_REF] Cootes | Building and using flexible models incorporating grey-level information[END_REF] a given model. At the basis of this activity is the growing availability of large databases of subjects and patients including biomedical signals and images as well as genetic information.

There is an additional dimension associated with a given model : the scale at which the anatomical, physical or physiological structure is described. With the development of new imaging modalities, it is now possible to observe the shape or function of most structures at the macroscopic (tissue), microscopic (cellular) levels and even in such cases to reveal the metabolic activity a the nanoscopic (molecular) scale.

There is a distinction between generic and specific models.

Generic models are useful for teaching [START_REF] Ryan | A virtual reality electrocardiography teaching tool[END_REF] and for training [START_REF] Delingette | Soft tissue modeling for surgery simulation[END_REF] medical residents, especially when they include a variability study between normals and pathological subjects. There is also a growing need for patient-specific models whose parameters are adjusted to the actual geometrical, physical or physiological [START_REF] Moreau-Villéger | Building maps of local apparent conductivity of the epicardium with a 2D electrophysiological model of the heart[END_REF] properties of the considered patient. This is motivated by the evolution of medical practice toward more quantitative and personalized decision processes for prevention, diagnosis and therapy.

In the next sections, following our proposed hierarchy (geometrical, physical, physiological), we describe a number of computational models used for image-guided, robot-assisted and simulated medical interventions before proposing some perspectives and challenges for future research.

II. GEOMETRICAL MODELING

A. Extracting anatomical information from images

Medical images are the primary source of information for the reconstruction of the human anatomy and are particularly well suited for the creation of patient-specific models. They can be generated by various imaging systems (MR, CT, US) and may have different dimensions (2D, 3D or 3D+t). The extraction of the geometrical description from those images requires the application of image segmentation algorithms that can generate a large variety of geometrical description of the targeted structure. The nature of the geometrical representation largely depends on the application. For instance, the volume rendering of a patient anatomy only requires to label the voxel belonging to the region of interest and to define its opacity function. For image-guided therapy a surface mesh of the structure is often sufficient while other therapy planning applications may need the definition of a volumetric mesh based on tetrahedra or hexahedra [START_REF] Zhang | Adaptive and quality quadrilateral/hexahedral meshing from volumetric data[END_REF].

There exist many pitfalls when creating meshes from medical images. For instance, the transformation of a segmented image into a surface mesh often requires using the Marching Cube algorithm [START_REF] Lorensen | Marching cubes: a high resolution 3D surface construction algorithm[END_REF] which tends to create a large number of poorly shaped triangles. Mesh decimation [START_REF] Schroeder | Decimation of triangles meshes[END_REF] is then applied to drastically reduce the number of elements and improve their shapes. Other meshing techniques [START_REF] Boissonnat | Provably good surface sampling and approximation[END_REF] are being developed that aim at directly producing a small number of well-shaped triangles from isosurfaces.

For volumetric meshes, the quality of the elements is especially of great importance since it influences the time of convergence and even the convergence of the simulation.

Even if hexahedral meshes usually lead to more accurate results than tetrahedral meshes [START_REF] Benzley | Comparison of all-hexahedral and all-tetrahedral finite element meshes for elastic and elasto-plastic analysis[END_REF] for the same number of nodes, their grid topology makes them ill-suited to represent complex shapes, unless a gross approximation of the surface is tolerated. The construction of patient-specific volumetric meshes can be produced from images in two ways. The first approach is bottom-up and consists in automatically filling a triangulated surface with tetrahedra [START_REF] Frey | Mesh Generation[END_REF], [START_REF] Alliez | Variational tetrahedral meshing[END_REF] or generating hexahedra from a set of voxels. The second approach is top-down and relies on a registration technique to deform a generic tetrahedral [START_REF] Sermesant | Deformable biomechanical models: Application to 4D cardiac image analysis[END_REF] or hexahedral mesh [START_REF] Chabanas | A 3D finite element model of the face for simulation in plastic and maxillo-facial surgery[END_REF] towards the segmented region in the image. The deformation of a generic model serves to automatically label mesh parts corresponding to specific boundary conditions, or specific material properties [START_REF] Sermesant | Deformable biomechanical models: Application to 4D cardiac image analysis[END_REF].

B. Hepatic Surgery Planning

Many state-of-the-art medical applications in clinical use rely on static anatomical models. For instance, we have developed [START_REF] Soler | Fully automatic anatomical, pathological, and functional segmentation from CT scans for hepatic surgery[END_REF] an automated hepatic surgery planning system that helps surgeons to choose a surgical strategy. This system can be used on a laptop during the intervention to improve the control of skill.

To acquire a knowledge of the patient anatomy, radiologists use CT-scans, which are acquired 60 seconds after the intravenous injection of a contrast medium. Theses images allow the visualization of hepatic tumors and contrasted vessels as well as the liver. We have developed a method which provides the segmentation of the following abdominal structures from a high resolution CT-scan within 15 minutes : skin, bones, lungs, heart, aorta, spleen, kidneys, gallbladder, liver, portal vein and hepatic tumors greater than 5 mm diameter (see Figure 2).

Note that other hepatic surgery planning systems [START_REF] Reitinger | Tools for augmented reality based liver resection planning[END_REF], [START_REF] Bourquain | A software assistant for preoperative planning in living-related liver transplantation and oncologic liver surgery[END_REF] have also been proposed.

C. Augmented Reality

Anatomical models are important not only for planning a therapy but also for guiding the execution of that therapy. This requires the additional task of relating the pre-operative plan coordinate system to the physical space where the patient actually lies. When the therapy is robot-assisted, this geometric correspondence and the anatomical knowledge of surgical field can be used to optimize the port placement [START_REF] Coste-Manière | Optimal planning of robotically assisted heart surgery: Transfer precision in the operating room[END_REF] or to define safety zones [START_REF] Starkie | Advances in active constraints and their application to minimally invasive surgery[END_REF] where the tip of the robot should not move. Similarly, Augmented Reality provides a convenient and efficient steering of the gesture by overlaying the pre-operative information onto real intra-operative views of the patient [START_REF]Proc. of the Scientific Workshop on Medical Robotics, Navigation and Visualization[END_REF].

Augmenting the surgeon's view with this information may lead to a better understanding of the pathology and decision process during the intervention [START_REF] Marescaux | Augmented reality system to guide radio-frequency tumor ablation[END_REF]. By tracking also the surgical tools in real time, or by robotizing them [START_REF] Taylor | Medical robotics in computer-integrated surgery[END_REF], one can also guide surgeons more efficiently along the planned path toward the virtual target.

Most applications developed so far are dedicated to neurosurgery and orthopedic surgery, since bones offer very reliable landmarks to register the virtual patient to the real one [START_REF] Grimson | An automatic registration method for frameless stereotaxy, image guided surgery and enhanced reality visualization[END_REF], [START_REF] Jolesz | Integration of interventional MRI with computer-assisted surgery[END_REF], [START_REF] Edwards | Design and evaluation of a system for microscope-assisted guided interventions (magi)[END_REF]. More recently, techniques based on the specific geometry of the intra-operative acquisition device were proposed for the chest or the abdomen, for instance to superimpose fluoroscopic images onto video images of the patient [START_REF] Mitschke | Recovering X-ray projection geometry for 3D tomographic reconstruction: Use of integrated camera vs. external navigation system[END_REF].

Other examples are given by the use of semi-transparent mirrors to display ultrasound [START_REF] Chang | Guidance of retrobulbar injection with real time tomographic reflection[END_REF], [START_REF] Stetten | C-mode real time tomographic reflection for a matrix array ultrasound sonic flashlight[END_REF] or CT slices [START_REF] Fichtinger | Needle insertion in CT scanner with image overlay -cadaver studies[END_REF], in such a way that images appear naturally registered to the patient.

However, these augmented reality techniques only superimpose intra-operative images onto video images of the patient.

To display pre-operative information, one needs to register those images to the patient geometry. For instance, we have developed [START_REF] Soler | Virtual reality and augmented reality in digestive surgery[END_REF] an augmented reality system for needle in- sertion guidance in hepatic radio-frequency ablation. This system superimposes in real-time 3D reconstructions from CT acquisitions and a virtual model of the needle on real video views of the patient, thanks to the 3D/2D registration of radioopaque markers stuck on the patient skin (see Figure 3).

D. Registration : a central issue

A careful evaluation of the performances of the registration algorithm [START_REF] Jannin | Validation of medical image processing in image-guided therapy[END_REF] is essential for augmented reality since some pre-operative information is used as a substitute for hidden information in video images: a registration error may lead the surgeon to miss the real target, with potential collateral damages.

Registration is basically an optimization procedure which aims at finding the best transformation parameters. It can be evaluated in terms of robustness, precision and accuracy.

First, one can distinguish between gross errors (convergence to wrong local minima) and small errors around the exact transformation. Robustness is the ability to converge to the right minimum, and may be quantified by the size of the basin of attraction, or the probability of false positives. Small errors around the right transformation might have a systematics trend (a bias), in which case the variance around the mean (characterizing the precision or reproducibility) is much smaller than the variance w.r.t. the ground truth (characterizing the accuracy). This distinction is important as one most often misses the ground truth to evaluate the performances of registration algorithms. In that case, one is left with the quantification of the precision, which can be significantly smaller than the real accuracy [START_REF] Pennec | Feature-based registration of medical images: Estimation and validation of the pose accuracy[END_REF]. Accuracy is an important concern in image guided therapy, and is indeed difficult to evaluate [START_REF] Nicolau | A complete augmented reality guidance system for liver punctures : first clinical evaluation[END_REF]. In some cases, however, one can certify the quality of the registration result in real-time, which allows to dynamically adapt the surgical procedure in order to reach the maximal safety [START_REF]An accuracy certified augmented reality system for therapy guidance[END_REF].

In terms of routine use of guidance systems in clinics, robustness can be even more important than accuracy, as rare but large errors may lead to much more dramatic consequences than a slightly larger average error. As a perfect robustness is not achievable, fault detection is thus of the upper importance. However, this is a very difficult problem: thanks to his anatomical interpretation of the image, a human observer immediately sees when the registration grossly fails, whereas the algorithm only knows that it is in a local minimum of a cost functional! Thus, providing a friendly user interface for the clinician to supervise and possibly interact to correct the registration results is often the key for the use of image guided systems in clinical environments [START_REF] Nicolau | An augmented reality & virtuality interface for a puncture guidance system: Design and validation on an abdominal phantom[END_REF].

Last but not least, registration algorithms have to be realtime in the image guided surgery context. By real-time, one means that the computation time should be significantly shorter than the typical motion time (a few minutes for the brain shift, a few tenths of second for the beating heart), the limiting factor being often the acquisition of the intraoperative images (e.g. 2 minutes for open MR in brain-shift compensation). In some cases, one might find some theoretical approximations that drastically speed-up the algorithm [START_REF] Granger | Multi-scale EM-ICP: A fast and robust approach for surface registration[END_REF].

However, one relies in most cases on the parallelization of the algorithm in order to reach a computation time compatible with the clinical constraints [START_REF] Stefanescu | A grid service for the interactive use of a parallel non-rigid registration algorithm of medical images[END_REF], [START_REF] Sermesant | A parallel implementation of non-rigid registration using a volumetric biomechanical model[END_REF].

III. BIOMECHANICAL/PHYSICAL MODELING

Anatomical models only provide a static geometrical representation of the patient anatomy and do not take into account the deformation of soft tissue that may occur before or during therapy. To address this issue, it is necessary to add a biomechanical model that can estimate soft tissue deformations under the application of known forces or displacements. The additional complexity of modeling may be used to improve the pre-operative planning of a therapy [START_REF] Hawkes | Tissue deformation and shape models in image-guided interventions: a discussion paper[END_REF], [START_REF] Schnabel | Validation of non-rigid image registration using finite element methods: application to breast MR images[END_REF] or to provide advanced surgical gesture training as illustrated in the sections below. Note that there are also important research efforts involving physical but non-mechanical phenomena for instance for the planning of radiotherapy or radiofrequency ablation [START_REF] Jain | A three-dimensional finite element model of radiofrequency ablation with blood flow and its experimental validation[END_REF].

A. Brain-Shift modeling in neurosurgery

The brain shift that occurs during a neuro-surgical intervention is the main source of intra-operative localization inaccuracies of pathologies (cerebral tumors,. . .). Indeed, a neurosurgeon establishes the surgical plan based on a preoperative MR image : any non-rigid motion of the brain between the pre-operative and the intra-operative configuration may lead to an error in the localization of the target. This clinical problem has motivated the study of the biomechanical behavior of the brain [START_REF] Miller | Constitutive modelling of abdominal organs[END_REF]. To model the brain motion after opening the dura, a number of authors [START_REF] Miga | Model-updated image guidance: Initial clinical experience with gravity-induced brain deformation[END_REF], [START_REF] Ferrant | Serial registration of intraoperative MR images of the brain[END_REF] have made the hypothesis that the loss of cerebro-spinal fluid causes a pressure field along the gravity direction (Archimedes principle). Furthermore, anatomical constraints (falx cerebri, skull) of the deformation field can be enforced with a biomechanical model of the brain discretized as a tetrahedral mesh since the relevant anatomical information can be extracted from MR images and enforced on the mesh. The predictive power of some models has been shown for instance for methods extrapolating displacement fields [START_REF] Ferrant | Serial registration of intraoperative MR images of the brain[END_REF], [START_REF] Skrinjar | Model-driven brain shift compensation[END_REF] from the cortex surface. Also, partial validation of brain shift models may be carried out [START_REF] Clatz | Patient specific biomechanical model of the brain: Application to parkinson's disease procedure[END_REF] by comparing computed displacements with those observed from intra-operative MR images.

Conversely, it interesting to stress that when validated, those biomechanical models can be used to register preoperative and intra-operative MR images [START_REF] Ferrant | Serial registration of intraoperative MR images of the brain[END_REF], [START_REF] Clatz | Robust non-rigid registration to capture brain shift from intra-operative MRI[END_REF] of the same patient in an efficient (near real-time) and realistic way.

Thus modeling brain tissue deformation leads to solving both direct (or predictive) and inverse (or estimation) problems and we propose to describe this global interaction as coupling biomechanical models with medical images. Building a surgical simulator requires the real-time and realistic simulation of soft-tissue deformations. More precisely the required refresh rate for user interaction is 30 Hz for visual feedback and more than 500 Hz for force-feedback, although the latter constraint can be alleviated with the addition of a local haptic model [START_REF] Balaniuk | Using fast local modeling to buffer haptic data[END_REF], [START_REF] Forest | Surface contact and reaction force models for laparoscopic simulation[END_REF]. Authors have proposed several soft tissue models trying to combine realism with efficiency.

B. Surgery Simulation

Discrete models such as spring-mass models are simple to implement but generally fail to produce realistic deformations [START_REF] Bianchi | Simultaneous topology and sti ness identification for mass-spring models based on FEM reference deformations[END_REF].

Combining finite element modeling and continuum mechanics has also been a popular approach which has provided advanced soft tissue simulation. Linear elasticity provides a simple and realistic framework when soft tissue undergoes small deformations. By exploiting the principle of superposition, one can precompute [START_REF] Delingette | Soft tissue modeling for surgery simulation[END_REF], [START_REF] Cotin | Real-time elastic deformations of soft tissues for surgery simulation[END_REF], [START_REF] James | ArtDefo accurate real time deformable objects[END_REF] the displacement field produced by applying a set of surface forces for an anisotropic and inhomogeneous material. This very efficient scheme (leading to a 500 Hz refresh rate or more) is however not compatible with mesh topological changes that are involved when simulating tissue resection. Tensor-mass models [START_REF] Cotin | A hybrid elastic model allowing real-time cutting, deformations and force-feedback for surgery training and simulation[END_REF] provide realistic volumetric deformations while authorizing topological changes (cuttings) with a complexity similar to spring-mass models.

To handle large deformations, one must introduce non-linear elastic models with the drawback of larger computational cost.

To reach real-time computation, authors have proposed nonlinear tensor-mass models [START_REF] Delingette | Soft tissue modeling for surgery simulation[END_REF], [START_REF] Picinbono | Non-Linear Anisotropic Elasticity for Real-Time Surgery Simulation[END_REF], multigrid [START_REF] Debunne | Dynamic realtime deformations using space and time adaptive sampling[END_REF] and corotational approaches [START_REF] Hauth | Corotational simulation of deformable solids[END_REF] each having its advantages and limitations in terms of computational efficiency, topological adaptation and deformation realism. Similarly, meshless methods [START_REF] Mueller | Meshless deformations based on shape matching[END_REF] are an interesting alternative to finite element methods where the location and the number of nodes can be optimized.

An important issue related to those models is the estimation of material properties such as the Young Modulus for linear elastic material or the stress-strain relationship for hyperelastic ones. There are three different sources of rheological data : exvivo testing where a sample of a tissue is positioned inside a testing rig [START_REF] Sakuma | In vitro measurement of mechanical properties of liver tissue under compression and elongation using a new test piece holding method with surgical glue[END_REF], [START_REF] Ottensmeyer | The effects of testing environment on the viscoelastic properties of soft tissues[END_REF]; in-vivo testing where a specific force and position sensing device is introduced inside the abdomen to perform indentation [START_REF] Nava | Determination of the mechanical properties of soft human tissues through aspiration experiments[END_REF], [START_REF] Chanthasopeephan | Measuring forces in liver cutting: New equipment and experimental results[END_REF]; Image-based elastometry from ultrasound, Magnetic Resonance Elastometry [START_REF] Kerdok | Truth Cube: Establishing Physical Standards for Real Time Soft Tissue Simulation[END_REF], [START_REF] Manduca | Magnetic resonance elastography: Non-invasive mapping of tissue elasticity[END_REF] or CTscan imaging.

There is no consensus on which method is best suited to recover meaningful material parameters. For instance, the rich perfusion of the liver affects deeply its rheology (the liver receives one fifth of the total blood flow at any time) and therefore it is still an open question whether ex-vivo experiments can assess the property of living liver tissue, even when specific care is taken to prevent the swelling or drying of the tissue. Conversely, data obtained from in-vivo experiments should also be considered with caution because the response may be location-dependent (linked to specific boundary conditions or non-homogeneity of the material) and the influence of the loading tool caliper on the deformation may not be well understood. Furthermore, both the respiratory and circulatory motions may affect in-vivo data.

Another difficult task when implementing a surgical simulation includes the detection and computation of collisions between a deformable model and rigid bodies [START_REF] Forest | Surface contact and reaction force models for laparoscopic simulation[END_REF] or between several deformable bodies [START_REF] France | A layered model of a virtual human intestine for surgery simulation[END_REF] (see Figure 5).

IV. PHYSIOLOGICAL MODELING

Up to now, we showed that geometrical and physical/biomechanical modeling of the human body provide representations allowing geometric reasoning, navigation, and various forms of simulated interactions including deforming and cutting soft tissues with realistic visual and haptic feedback. However the tissues and organs still behave like passive material, and it is necessary to go one step further to model the active properties of living tissues and the dynamic nature of normal or pathological evolving processes. We illustrate this level of modeling with two examples related to the modeling of the electro-mechanical activity of the heart and the growth of brain tumors. Additional examples related to the modeling of physiological flows can be found in [START_REF] Quarteroni | Mathematical modeling and numerical simulation of the cardiovascular system[END_REF], [START_REF] Boissonnat | From arteriographies to computational flow in saccular aneurisms: the INRIA experience[END_REF], [START_REF] Ztonyi | Real-time synthesis of bleeding for virtual hysteroscopy[END_REF]. (from [START_REF] Delingette | Soft tissue modeling for surgery simulation[END_REF], [START_REF] Cotin | Real-time elastic deformations of soft tissues for surgery simulation[END_REF], [START_REF] Cotin | A hybrid elastic model allowing real-time cutting, deformations and force-feedback for surgery training and simulation[END_REF] and [START_REF] Forest | Removing tetrahedra from manifold tetrahedralisation : application to real-time surgical simulation[END_REF], [START_REF] Forest | Surface contact and reaction force models for laparoscopic simulation[END_REF]).

for medical purposes. The model reproduces the electrical depolarization and repolarization of the cardiac tissues through a set of macroscopic reaction-diffusion equations initially proposed by Fitzugh and Nagumo [START_REF] Fitzhugh | Impulses and physiological states in theoretical models of nerve membrane[END_REF] and further refined by Aliev and Panvilov [81]. This electrical activity, which can be synchronized with the actual ECG (electro-cardiogram) of the patient, creates a mechanical contraction followed by an active relaxation which are modeled by a set of partial differential equations initially proposed by Bestel, Clement and Sorine [START_REF] Bestel | A biomechanical model of muscle contraction[END_REF].

It was shown in [START_REF] Sermesant | An electromechanical model of the heart for image analysis and simulation[END_REF] that this model could be interactively adjusted to the actual geometrical, mechanical or electrical properties of patient's heart through the use of conventional or tagged MR images and some in vivo electrophysiological measurements. The average direction of the myocardium fibers is also integrated into this model, because it plays an important role in both the electrical and mechanical modeling. Even more interestingly, the model can then be used to study the effect of modifying locally some electrical or mechanical properties in order to better predict the effect of a therapy or the evolution of a pathology [START_REF] Sermesant | Cardiac function estimation from mri using a heart model and data assimilation: Advances and difficulties[END_REF], [START_REF] Sermesant | Simulation of the electromechanical activity of the heart using xmr interventional imaging[END_REF] (see Figure 6). We believe that this kind of models of a dynamic organ could be used in the future to better plan or train a number of medical skills for instance in radiofrequency ablation procedure, for the positioning of stimulation probes of pacemakers, and maybe in the long term future to simulate beating heart with the potential use of robotics [START_REF] Ginhoux | Beating heart tracking in robotic surgery using 500 hz visual servoing, model predictive control and an adaptive observer[END_REF]. The ejection fraction decreases from 65% to 55%.

B. Tumor growth

The second example is related to the modeling of the growth of brain tumors. A joint action between INRIA, a Nice Hospital (Centre Antoine Lacassagne) and the Brigham and Women's hospital (Boston) has led to the development of a three level model [START_REF] Clatz | In silico tumor growth: Application to glioblastomas[END_REF].

The first level includes the geometrical model of a patient's head, including the skull, the brain parenchyma (grey and white matter) and the Cerebro-Spinal Fluid (CSF). The shape of each region is acquired through a conventional MR exam.

In addition, the direction of the main white matter fibers is also acquired either through Diffusion Tensor Imaging (DTI), or using average directions provided by a brain atlas. A second level includes the modeling of the biomechanical properties of the brain tissues. Because we are considering small deformations only, this model is a linear elastic one, and implements the boundary conditions imposed by the bony (skull) and fluid (ventricles) structures surrounding the brain parenchyma. The third level is an evolution model of the tumoral tissues, which is based on a set of reaction-diffusion equations describing the proliferation and diffusion of malignant cells [START_REF] Swanson | Virtual brain tumours (gliomas) enhance the reality of medical imaging and highlight inadequacies of current therapy[END_REF].

An original point is the coupling of the third level with the previous two levels in order to create a realistic deformation of the brain tissues (also called mass effect) produced by the tumor growth. The direction of the white matter fibers plays an important role in the highly anisotropic diffusion of cancerous cells (see Figure 7). We believe that this type of model will prove to be quite useful to predict the evolution of a pathology from a number of observations made at a limited number of time points. It will help the guidance of medical strategies (e.g. radiotherapy, surgery) in predicting more precisely the boundaries of • Real-time Coupling of Models with Observations. This requires very fast inversion of realistic computational models. This will improve robot-assisted image-guided therapy in the presence of patient motion (brain shift, respiration, cardiac motion, etc). Figure 8 previews the combination of augmented reality (geometrical models) with medical robotics.

• Miniaturized Robotics The advances in robotics research, in particular the design of new generations of highly accurate, easy-to-use miniature robotic guidance systems (e.g. the spine-assist system [START_REF] Shoham | Bone-mounted miniature robot for surgical procedures: Concept and clinical applications[END_REF]) might also contribute to the successful combination of computational models and medical imaging.

• Microscopic Imaging. New in vivo microendoscopy tech- niques [START_REF] Macaulay | In vivo pathology: microendoscopy as a new endoscopic imaging modality[END_REF], [START_REF] Goualher | Towards optical biopsies with an integrated fibered confocal fluorescence microscope[END_REF] providing structural information on the tissues at the cellular level should also open, combined with robotics and computational models, new avenues for improving diagnosis and therapy.

Fig. 2 .

 2 Fig.2. Reconstructed liver from the patient CT scan image (from[START_REF] Soler | Fully automatic anatomical, pathological, and functional segmentation from CT scans for hepatic surgery[END_REF]).

Fig. 3 .

 3 Fig. 3. (Left) Augmented reality based on 3D radio-opaque markers; (Right) Augmented view of patient at the tip of the needle (from [43]).

  Surgery simulation aims at reproducing the visual and haptic senses experienced by a surgeon during a surgical procedure, through the use of computer and robotics systems. The medical scope of this technology is linked with the development of minimally invasive techniques especially videoscopic surgery (endoscopy, laparoscopy,...) and possibly telesurgery. Because of the limitations of current training methods, there is a large interest in developing video-surgery simulation software for providing efficient and quantitative gesture training systems. Indeed, such systems should bring a greater flexibility in the training by providing scenarios including different types of pathologies. Furthermore, by creating patient-specific models, surgery simulation allows surgeons to verify, optimize and rehearse the surgical strategy of a procedure on a given patient.

Fig. 5 .

 5 Fig. 5. Simulation of an intervention on a human intestine [78] that involves complex handling of self-collisions (courtesy of M-P. Cani).

A

  . Cardiac modeling During the past 5 years, a scientific INRIA consortium regrouping internal and external research teams 1 has developed an electro-mechanical model of the cardiac ventricles

Fig. 4 .

 4 Fig. 4. (Left) A force feedback system suited for surgery simulation; (Right) View of the simulated hepatic resection involving linear-elastic materials

Fig. 6 .

 6 Fig. 6. Predicting the effect of pathologies with an electro-mechanical model of the heart [8]; (Left) Definition of a basal left epicardial infarcted zone (no conduction and no contraction in that area) shown in blue; (Right) Curves of ejected volume during a cardiac cycle without (dashed) and with (solid) infarct.

Fig. 7 .

 7 Fig. 7. (Top) Overview of the evolution model of cancerous cells which takes into account the anisotropic diffusion along white matter fibers; (Bottom) Result of the tumor growth simulation on a brain slice; (a) Initial MR T2 image of the patient with lines of constant tumor density; (b) View of the corresponding MR T2 slice (after rigid registration) six months later; (c) Lines of constant tumor density predicted by the tumor growth model [86] after 6 months of evolution.

Fig. 8 .

 8 Fig.8. By combining Augmented reality (left) and robotics (right) accuracy of image guided surgery may be improved and this may lead to additional levels of surgical automation (from[START_REF] Marescaux | Augmented reality system to guide radio-frequency tumor ablation[END_REF]).

cf. CardioSense3D URL http://www.inria.fr/CardioSense3D/


X Pennec

P Cachier
email: pascal.cachier@sophia.inria.fr

N Ayache

Tracking brain deformations in time-sequences of 3D US images

Keywords: Non rigid registration, 3D ultrasound images, Tracking

During a neurosurgical intervention, the brain tissues shift and warp. In order to keep an accurate positioning of the surgical instruments, one has to estimate this deformation from intra-operative images. We present in this article a feasibility study of a tracking tool based on intra-operative 3D ultrasound (US) image sequences. The automatic processing of this kind of images is of great interest for the development of innovative and low-cost image guided surgery tools. The difficulty relies both in the complex nature of the ultrasound image, and in the amount of data to be treated as fast as possible.

Introduction

The use of stereotactic systems is now a quite standard procedure for neurosurgery. However, these systems are only accurate under the assumption that the skull and the brain move together as a unique rigid body during surgery. In practice, relative motion of the brain with respect to the skull (also called brain shift) occurs, mainly due to tumor resection, cerebrospinal fluid drainage, hemorrhage or even the use of diuretics. Furthermore, this motion is likely to increase with the size of the skull opening and the duration of the operation.

Over the last years, the development of real-time 3D ultrasound (US) imaging has revealed a number of potential applications in image-guided surgery as an alternative approach to open MR and intra-interventional CT. The major advantages of 3D US over existing intra-operative imaging techniques are its comparatively low cost and simplicity of use. However, the automatic processing of US images has not gained the same degree of development as other medical imaging modalities, probably due to the low signal-to-noise ratio of US images.

Context

We present in this article a feasibility study of a tracking tool for brain deformations based on intra-operative 3D ultrasound (US) image sequences. This work was performed within the framework of the European project ROBO-SCOPE, a collaboration between The Fraunhofer Institute (Germany), Fokker Control System (Netherlands), Imperial College (UK), INRIA (France), ISM-Salzburg and Kretz Technik (Austria). The goal of the whole project is to assist neuro-surgical operations using real-time 3D ultrasound images and a robotic manipulator arm (fig. 1). The operation is planned on a pre-operative MRI (MR1) and 3D US images are acquired during surgery to track in real time the deformation of anatomical structures. The first US image (US1) is acquired with dura mater still closed and a rigid registration with the preoperative MR is performed. This allows to relate the MR and the US coordinate systems and possibly to correct for the distortions of the US acquisition device. Then, brain deformations are tracked in the time-sequence of per-operative US images. From these deformations, one can update the preoperative plan and synthetize a virtual MR image that matches the current brain anatomy.

MR/US registration

The idea of MR/US registration is already present in (Trobaugh et al., 1994a;[START_REF] Trobaugh | Three-dimensional imaging with stereotactic ultrasonography[END_REF] where the US probe is calibrated (i.e. registered to the surgical space) and then tracked using an optical device. The MR/US registration is then obtained using the registration of the MR image to the surgical space by standard stereotactic neurosurgical procedures. [START_REF] Richard | A low-cost PCIbus-based ultrasound system for use in image-guided neurosurgery[END_REF] improved this method by designing a real-time low-cost US imaging system based on a PCI bus. In [START_REF] Pagoulatos | Interactive 3-D Registration of Ultrasound and Magnetic Resonance Images Based on a Magnetic Position Sensor[END_REF], the tracking of the US probe is performed with a DC magnetic position sensor. In [START_REF] Erbe | 3D-Ultrasonography and Image Matching for Detection of Brain Shift During Intracranial Surgery[END_REF], the registration is performed by interactively delineating corresponding surfaces in all images and a visual rigid fitting of the surfaces using a 6D space-mouse. In [START_REF] Hata | Registration of Ultrasound echography for Intraoperative Use: A Newly Developed Multiproperty Method[END_REF], the outlines of the 2D US image is registered to the MR surface using a Chamfer matching technique. All these techniques only perform a rigid registration of the MR and the US images.

For a non rigid registration (i.e. a brain shift estimation), we have to turn to [START_REF] Gobbi | Ultrasound probe tracking for real-time ultrasound/MRI overlay and visualization of brain shift[END_REF][START_REF] Gobbi | Ultrasound/MRI overlay with image warping for neurosurgery[END_REF][START_REF] Comeau | Intraoperative ultrasound for guidance and tissue shift correction in image-guided neurosurgery[END_REF], where the 2D US probe is still optically and rigidly tracked but the corresponding MR slice is displayed to the user who marks corresponding points on MR and US slices. Then, a thin plate spline warp is computed to determine the brain shift. This method is also developed in [START_REF] Bucholz | The correction of stereotactic inaccuracy caused by brain shift using an intraoperative ultrasound device[END_REF] with the possibility of using 3D US images and a deformation computed using a spring model instead of splines. More recently, Ionescu et al [START_REF] Ionescu | Automated Registration of Ultrasound with CT Images: Application to Computer Assisted Prostate Radiotherapy and Orthopedics[END_REF] registered US with Computed Tomography (CT) data after automatically extracting contours from the US using watershed segmentation. In these studies, there is no processing of a full time sequence of US images : the brain shift estimation is limited to a few samples at given time-points as the user interaction is required at least to define the landmarks.

Recently, an automatic rigid registration of MR and US images was presented [START_REF] Roche | Rigid regis-tration of 3D ultrasound with MR images: a new approach combining intensity and gradient information[END_REF][START_REF] Roche | Generalized Correlation Ratio for Rigid Registration of 3D Ultrasound with MR Images[END_REF]Pennec et al., 2001a). This work is based on image intensities and does not rely on feature extraction. However, the estimated motion remains limited to rigid or possibly affine transformations. Up to our knowledge, only [START_REF] King | Baysian estimation of intra-operative deformation for imageguided surgery using 3-d ultrasound[END_REF] deals with an automatic nonrigid MR/US registration: the idea is to register a surface extracted from the MR image to the 3D US image using a combination of the US intensity and the norm of its gradient in a Bayesian framework. The registration is quite fast (about 5mn), even if the compounding of the 3D US and the computation of its gradient takes about one hour. However, experiments are presented only on phantom data and our experience (see section 3) is that real US images may lead to quite different results.

Tracking methods in sequences of US images

Since non-rigid MR/US registration is a difficult problem, we chose to split it into two subproblems: first a rigid MR/US registration is performed with dura matter still closed (there is no brain shift yet), for instance using the approach of [START_REF] Roche | Generalized Correlation Ratio for Rigid Registration of 3D Ultrasound with MR Images[END_REF]. Then we look for the non-rigid motion within the US time-sequence. In the literature, we found a small number of articles on the registration of 3D US images. [START_REF] Strintzis | Maximum Likelihood Motion Estimation in Ultrasound Image Sequences[END_REF]) use a maximum-likelihood approach to deduce a similarity measure for ultrasound images corrupted by Rayleigh noise and a block-matching strategy to recover the rigid motion. In [START_REF] Rohling | Automatic registration of 3-D ultrasound images[END_REF], the correlation of the norm of the image gradient is used as the similarity measure to rigidly register two US images in replacement of the landmark-based RANSAC registration of [START_REF] Rohling | Three-Dimensional Spatial Compounding of Ultrasound Images[END_REF]. However, these methods only deal with rigid motion and consider only two images, eluding the tracking problem.

One has to move to cardiac application to find some real tracking of nonrigid motion in US images. In [START_REF] Papademetris | 3D cardiac deformation from ultrasound images[END_REF], the endo-and epicardial surfaces are interactively segmented on each 2D image plane. Then, a shape-memory deformable model determines the correspondences between the points of the 3D surfaces of successive images. These correspondences are used to update an anisotropic linear elastic model (finite element mesh). The approach is appealing but relies once again on an interactive segmentation. In [START_REF] Sanchez-Ortiz | Automatic 3D echocardiographic image analysis[END_REF], a combination of feature point extraction (phasebased boundaries) and a multi-scale fuzzy clustering algorithm (classifying the very low intensities of intra-ventricular pixels) is used to segment the surface of the left ventricular cavity. This process is done in 2D+T and then reconstructed in 3D. Thus it exploits the whole sequence before tracking the motion itself, which is not possible for our application. These two methods are well suited for the shape of the cardiac ventricle using dedicated surface models. If they could be adapted to the brain ventricles, it seems difficult to extend them to the tracking of the volumetric deformations of the whole brain.

Intensity based non-rigid registration algorithms

Since feature or surface extraction is especially difficult in US images, we believe that an intensity-based method can more easily yield an automatic algorithm. Over recent years, several non-rigid registration techniques have been proposed. [START_REF] Bajcsy | Multiresolution Elastic Matching[END_REF]) differentiated the linear correlation criterion and used a fixed fraction of its gradient as an external force to interact with a linear elasticity model. [START_REF] Christensen | Volumetric Transformation of Brain Anatomy[END_REF] show that the linear elasticity, valid for small displacements, cannot guarantee the conservation of the topology of the objects as the displacements become larger: the Jacobian of the transformation can become negative. Thus, he proposed a viscous fluid model of transformations as it can handle larger displacement. This model is also linearized in practice. [START_REF] Bro-Nielsen | Medical image registration and surgery simulation[END_REF] started from the fluid model of Christensen and used the linearity of partial derivative equations to establish a regularization filter, several order of magnitude faster than the previous finite element method. He also justified his forces as the differential of the sum of square intensity differences criterion, but he still used a fixed fraction of this gradient, and shows that Gaussian smoothing is an approximation of the linear elastic model. Some authors [START_REF] Maintz | General Multimodal Elastic Registration based on Mutual Information[END_REF] tried to apply to non-rigid registration some criteria developed for rigid or affine matching using block-matching techniques. However, these criteria require a minimal window size, thus limiting the resolution of the result. Moreover, the regularization of the displacement field is usually implicit, i.e. only due to the integration of the criterion over the window, which means that it is difficult to explicitly control the regularity of the sought transformation.

Recently, [START_REF] Thirion | Image matching as a diffusion process: an analogy with Maxwell's demons[END_REF] proposed to consider non rigid registration as a diffusion process. He introduced in the images entities (demons) that push according to local characteristics of the images in a similar way Maxwell did for solving the Gibbs paradox in thermodynamics. The forces he proposed were inspired from the optical flow equations. This algorithm is increasingly used in several teams as reported by [START_REF] Dawant | Brain atlas deformation in the presence of large space-occupying tumors[END_REF][START_REF] Bricault | Registration of Real and CT-Derived Virtual Bronchoscopic Imag es to Assist Transbronchial Biopsy[END_REF][START_REF] Webb | Automatic Detection of Hippocampal Atrophy on Magnetic Resonnance Images[END_REF][START_REF] Prima | Automatic Analysis of Normal Brain Dissymmetry of Males and Females in MR Images[END_REF]. In [START_REF] Pennec | Understanding the "demon's algorithm": 3D non-rigid registration by gradient descent[END_REF][START_REF] Cachier | Fast non-rigid matching by gradient descent: Study and improvements of the "demons" algorithm[END_REF], we investigated the non-rigid registration using gradient descent techniques. Differentiating the sum of square intensity differences criterion (SSD), we showed that the demons forces are an approximation of a second order gradient descent on this criterion. The same gradient descent techniques were applied to a more complex similarity measure in [START_REF] Cachier | 3D non-rigid registration by gradient descent on a gaussian-windowed similarity measure using convolutions[END_REF]: the sum of Gaussian-windowed local correlation coefficients (LCC).

Overview of the article organization

In this article, we develop an automatic intensity-based non-rigid tracking algorithm suited for real-time US images sequences, based on encouraging preliminary results reported in [START_REF] Pennec | Understanding the "demon's algorithm": 3D non-rigid registration by gradient descent[END_REF][START_REF] Cachier | Fast non-rigid matching by gradient descent: Study and improvements of the "demons" algorithm[END_REF][START_REF] Pennec | Tracking brain deformations in time sequences of 3D US images[END_REF]Pennec et al., 2001a). We first present the registration method for two US images. We detail in section 2.1 a new parameterized deformation field. Then, we define in sections 2.2 and 2.3 the similarity and regularization energies, which are optimized in section 2.4 using a gradient descent algorithm. We show in section 2.6 how to turn the registration of two images into a tracking algorithm. In section 3, we present qualitative results of the tracking algorithm on a sequence of US images of a phantom, and quantitative results on a small sequence of US images of a dead pig brain with a simulated cyst. Our results tend to justify the choices of the similarity energy and of the model of deformations, particularly with the longer term goal of achieving a real-time tracking system.

The tracking algorithm

When analyzing the problem of tracking the brain deformation in 3D US timesequences, we made the following observations. First, deformations are small between successive images in a real-time sequence, but they are possibly large deformations around the surgical tools with respect to the pre-operative image. Thus, the transformation space should allow large deformations, but only small deformations have to be retrieved between successive images. Second, there is a poor signal to noise ratio in US images and the absence of information in some areas. However, the speckle (inducing localized high intensities) is usually persistent in time and may produce reliable landmarks for successive images [START_REF] Meunier | Ultrasonic Texture Motion Analysis: Theory and Simulation[END_REF]. As a consequence, the transformation space should be able to interpolate in areas with few information while relying on high intensity voxels for successive images registration. Last but not least, the algorithm is designed in view of a real-time registration during surgery, which means that, at equal performances, one should prefer the fastest method.

Following the encouraging results obtained in [START_REF] Pennec | Understanding the "demon's algorithm": 3D non-rigid registration by gradient descent[END_REF][START_REF] Cachier | Fast non-rigid matching by gradient descent: Study and improvements of the "demons" algorithm[END_REF] for the intensity based non-rigid registration of two 3D US images, we adapt in this section the method in four different directions, according to the previous observations. We first look for more robust free-form transformations. Then we compare different image similarity criteria and different optimization strategies. Finally, we transform the registration algorithm into a tracking tool suited for time sequences.

Parameterization of the transformation

The brain shift is only a small deformation, but the introduction of surgical tools and the removal of tissues near the region of interest may locally introduce some large deformations. Simple transformations, like rigid or affine ones, can be represented by a small number of parameters (resp. 6 and 12 in 3D). When it comes to free-form deformations, we need to specify the coor-dinates T (x) of each point x of the image after the transformation. Such a non-parametric transformation is usually represented by its displacement field U (x) = T (x) -x (or U = T -Id), sampled at each voxel. This strategy proved to be successful in textured enough regions but induces convergence problems in large uniform areas (as it is the case in the phantom sequence of section 3.1) because the propagation of regularization constraints is very slow.

We found that a re-parameterization of the transformation was necessary to promote a better conditioning of the problem. We previously had a displacement t i for each voxel position x i . Now, t i represents a parameter of a smooth transformation defined by:

T (t 1 , ...t n )(x) = i t i .G σ (x -x i ) (1)
Note that when σ goes to 0, the parameterization tends toward the previous parameterization. Moreover, this parameterization can still interpolate any value at each site x i . The transformation being described as being a sum of Gaussians, rather than a sum of Diracs, the gradient descent algorithm uses the derivatives of the similarity with respect to the displacement of an entire group of voxels, which is more robust to noise, and will lead to a faster propagation of regularity constraints in uniform intensity areas.

In this article, we used a site x i at each voxel of the image. One could think of reducing the number of sites to decrease the number of parameters and go toward smoother transformations. However, we observed that this is not equivalent to the smoothing performed in section 2.3: the Gaussian introduced here parameterizes the width of the neighborhood around voxel x i for which the voxel intensities will have an influence on the transformation parameter t i at x i . Thus, it can be seen as a regularization of the similarity energy landscape, as we will see in Eq. 6, and not as a regularization of the transformation. Thus, reducing the number of sites (at a fixed Gaussian width) would only reduce the resolution of the transformation.

Similarity energy

Even if there is a poor signal to noise ratio in US images, the speckle is usually persistent in time and may produce reliable landmarks within the time-sequence [START_REF] Meunier | Ultrasonic Texture Motion Analysis: Theory and Simulation[END_REF]. Hence, it is desirable to use a similarity measure which favors the correspondence of similar high intensities for the registration of successive images in the time-sequence. First experiments presented in [START_REF] Pennec | Understanding the "demon's algorithm": 3D non-rigid registration by gradient descent[END_REF][START_REF] Cachier | Fast non-rigid matching by gradient descent: Study and improvements of the "demons" algorithm[END_REF] indicated that the simplest one, the sum of square differences (SSD), could be suited. Let I be the reference image and J • T the transformed image to register; the criterion to minimize is:

SSD(T ) = (I -J • T ) 2
In [START_REF] Cachier | 3D non-rigid registration by gradient descent on a gaussian-windowed similarity measure using convolutions[END_REF], we developed a more complex similarity measure: the sum of Gaussian-windowed local correlation coefficients (LCC). Let G f denote the convolution of f by the Gaussian function G. We define the local mean by Ī = (G I), the local variance by σ 2 I = G (I -Ī) 2 , and the local correlation of I and

J • T by LC(T ) = G (I -Ī)(J • T -J • T ) .
Then, the global criterion to maximize is the sum of the local correlation coefficients:

LCC(T ) = LC(T ) σ I .σ J•T
We showed that this criterion can be differentiated up to the second order using only recursive Gaussian convolutions which are very fast and in a time independent of the standard deviation of the Gaussian. Thus it may be optimized using a gradient descent like the SSD criterion.

We run most of the experiments presented in section 3 with the LCC criterion and we did not find significant differences with the results of the algorithm using the SSD. However, the LCC is still around 2 times slower than the SSD. Since the computation time of the US-US non-rigid registration is a key issue for real-time motion tracking, we preferred to keep the SSD criterion. We believe that this choice is justified anyway for the registration of successive images in the time sequence, but it could be reconsidered for the update of the global deformation (transformation from the first image to the current one, see section 2.6) if the sequence was to present some important intensity changes along time.

Regularization energy

In non-rigid registration, there is a trade-off to find between the similarity energy, reflected by the visual quality of the registration, and the smoothing energy, reflected by the regularity of the transformation (the term "regularity" should be taken in its broadest sense, since the smoothing energy may allows occasional discontinuities in the displacement field [START_REF] Hellier | Medical Image Registration with Robust Multigrid Techniques[END_REF]).

In the regularization theory framework, one minimizes the weighted sum of the energies: E sim + λ.E reg . This formulation has proven to be successful for data approximation, and has been used for various approaches of non-rigid registration algorithms [START_REF] Ferrant | 3D Image Matching using a Finite Element Based Elastic Deformation Model[END_REF]. However, there is an important difference between data approximation and image registration. In data approximation, both energies measure different properties of the same object (the similarity and the smoothness of the data), while the two energies relate to different objects in image registration (the intensities of the images for the matching energy and the transformation for the regularization energy). Thus, one has to find a non linear tradeoff between the two energies.

Another widely spread method attempts to separate the image measure from the transformation measure, and could be compared with the approach of game theory. It consists in alternatively decreasing the similarity energy and the smoothing energy. This approach is chosen in many block-matching algorithms [START_REF] Ourselin | Block matching: a general framework to improve robustness of rigid registration of medical images[END_REF] and in some optical-flow-based techniques [START_REF] Thirion | Image matching as a diffusion process: an analogy with Maxwell's demons[END_REF]. In view of a real-time system, this is particularly well suited for the stretch energy (or membrane model) E reg = ||∇T || 2 = Tr(∇T.∇T T ) as the associated Euler-Lagrange evolution equation corresponds to the heat propagation in a homogeneous material. Thus, one step of gradient descent corresponds to convolution of the transformation by a Gaussian with a standard deviation linked to the time step of the gradient descent [START_REF] Morel | Variational Methods in Image Segmentation[END_REF]. This way, we get a simple regularization by a Gaussian smoothing of the transformation parameters t i with a smoothing parameter (the σ T of this Gaussian) that has a physical meaning.

In summary, the algorithm will alternatively perform one step of gradient descent on the similarity energy E sim and one step of transformation smoothing by Gaussian filtering of standard deviation σ T .

Minimizing the similarity energy for a free-form deformation

Let T be the current estimation of the transformation and (∇ J • T )(x) (resp. (H J • T )(x)) be the transformed gradient (resp. Hessian) of the image J. A perturbation by a displacement field u(x) gives the following Taylor expansion:

(J • (T + u))(x) = (J • T )(x) + (∇ J • T ) T .u(x) + 1 2 u(x) T .(H J • T ).u(x)
Thus, the Taylor expansion of the criterion is:

SSD(T + u) = SSD(T ) + 2 (J • T -I) .(∇ J • T ) T .u + ((∇ J • T ) T .u) 2 + (J • T -I) .u T .(H J • T ).u + O(||u|| 2 )
where ||u|| 2 = x u(x) 2 .dx is the L 2 norm of the small perturbation. As, by definition, x f (x) T .u(x).dx is the dot product of f and u in the space of square-integrable functions, we get by identification:

∇ SSD (T ) = 2(J • T -I).(∇ J • T ) (2) H SSD (T ) = 2(∇ J • T ).(∇ J • T ) T + 2(J • T -I).(H J • T ) (3)
Let us now approximate the criterion by its tangential quadratic form at the current transformation T . We get the following first order approximation of the criterion gradient: ∇ SSD (T + u) ∇ SSD (T ) + H SSD (T ).u

Assuming that the Hessian matrix of the criterion is positive definite, the minimum is obtained for a null gradient, i.e.: u = -H (-1) SSD (T ).∇ SSD (T ). This formula require to invert the Hessian matrix H SSD (T ) at each point x of the image. To speed up the process, we approximate this matrix by the closest scalar matrix (for the L 2 norm on the matrix vector space):

H SSD (T ) Tr (H SSD (T )) n .Id = ∇ J • T 2 + (J • T -I).(∆ J • T ) 3 .Id
where n is the space dimension (3 for us). Using this approximation, we get the following adjustment vector field:

u -3.(J • T -I).(∇ J • T ) ||∇ J • T || 2 + (J • T -I).(∆ J • T ) (4)
In fact, when minimizing the reverse SSD criterion (I • T (-1) -J) 2 , one finds that the optimal adjustment is given by [START_REF] Pennec | Understanding the "demon's algorithm": 3D non-rigid registration by gradient descent[END_REF][START_REF] Cachier | Fast non-rigid matching by gradient descent: Study and improvements of the "demons" algorithm[END_REF]:

T = T • (Id + u ) with u = 3.(I -J • T ).∇ I ||∇ I || 2 + (I -J • T ).∆ I
which justifies the empirical force used by Thirion's demons:

v = (I -J • T ).∇ I ||∇ I || 2 + α.(I -J • T ) 2
In practice, we have modified the Newton optimization scheme described above into a Levenberg-Marquardt method where the adjustment vector field is given at each step by u = -(λ.Id+H SSD ) (-1) .∇ SSD . Dropping the second order terms in the Hessian, we are left with:

u = -3.(J • T -I)/(||∇ J • T || 2 + λ 2 ).(∇ J • T ) (5)
The parameter λ performs a tradeoff between a first order gradient descent (λ 1 means that we don't trust the Hessian matrix and we simply go along the gradient with a small time-step) and a second order gradient descent (λ 1 means that we use our simplified Hessian matrix). At each step, λ is divided by a fixed value α (typically 5) if the similarity criterion decreased, and the criterion is re-estimated with λ multiplied by α otherwise until the criterion decreases.

Minimizing the similarity energy for the new type of transformations

We now detail the differences induced by our new parameterization of the freeform transformation on the SSD criterion. Using the Gaussian parameterization of the transformation (eq. 1), t i is now a parameter of the transformation. Let G (x i ,σ) f denote the convolution by a Gaussian of variance σ centered at x i . Deriving the SSD w.r.t. this parameter gives:

∇ SSD (T ) = 2 G (x i ,σ) (J • T -I).(∇J • T ) (6) H SSD (T ) = 2 G 2 (x i ,σ) (∇J • T ).(∇J • T ) T + (J • T -I).(HJ • T ) (7)
Thus, the Gaussian parameterization acts as a smoothing on the gradient and Hessian of the energy. Therefore, it will be more robust and may escape from previous local minima. The minimization is performed as above with a Levenberg-Marquardt method using these regularized version of the energy derivatives.

From the registration to the tracking algorithm

In the previous sections, we studied how to register two US images together. We now have to estimate the deformation of the brain between the first image (since the dura mater is still closed, it is assumed to correspond to the preoperative brain) and the current image of the sequence. One could think of registering directly U S 1 (taken at time t 1 ) and U S n (at time t n ) but the deformations could be quite large and the intensity changes important. To constrain the problem, we need to exploit the temporal continuity of the deformation.

First, assuming that we already have the deformation T U S (n) from image U S 1 to U S n , we register U S n with the current image U S n+1 , obtaining the transformation dT U S (n). If the time step between two images is short with respect to the deformation rate (which should be the case in real-time sequences at a rate ranging from 1 to 5 images per second), this registration should be easy.

Moreover, the intensity changes should be small. For this step, we believe that the SSD criterion is well adapted.

Then, composing with the previous deformation, we obtain a first estimation of T U S (n + 1) dT U S (n) • T U S (n). However, the composition of deformation fields involves interpolations and just keeping this estimation would finally lead to a disastrous cumulation of interpolation errors:

T U S (n + 1) = dT U S (n) • dT U S (n -1) . . . dT U S (2) • dT U S (1)
Moreover, a small systematic error in the computation of dT U S (n) leads to a huge drift in T U S (n) as we go along the sequence. Thus, we only use dT U S (n) • T U S (n) as an initialization for the registration of U S 1 to U S n . Starting from this position, the residual deformation should be small (it corresponds to the correction of interpolation and systematic error effects) but the difference between homologous point intensities might remain important. In this case, the LCC criterion might be better than the SSD one despite its worse computational efficiency.

One of the main consequences is that the first US image will have to be of very high quality since it will be the only reference for tracking deformations along the whole sequence. One possibility consists in acquiring several images of the still brain in order to compute a mean image of better quality. Another possibility consists in performing some anisotropic diffusion on U S 1 to improve its quality.

Experiments

In this section, we present qualitative results of the tracking algorithm on a sequence of US images of a phantom, and quantitative results on a small sequence of US images of a dead pig brain with a simulated cyst. Experiments were performed using the SSD and the LCC criterion without significant differences in the results. Since the LCC is around 2 times slower than the SSD, we present here results and computation times for the SSD criterion. All 3D-US images were acquired using a commercial 3D-US volume scanner Voluson 530 D from Kretz Technology (4-9 MHz, 90 degrees aperture).

A Phantom study

Within the ROBOSCOPE project, an MR and US compatible phantom was developed by Prof. Auer and his colleagues at ISM (Austria) to simulate brain deformations. It is made of two balloons, one ellipsoid and one ellipsoid with a "nose", that can be inflated with known volumes. Each acquisition consists in one 3D MR image and one 3D US image (see Fig. 3 for an example). The goal is to use the US sequence to track the deformations and compute the corresponding virtual MR images from the first MR image. Then, the original MR images can be used to assess the quality of the tracking.

Since the US probe cannot enter the MR machine, it was removed for the MR acquisitions. Thus, we had to compensate for the apparent motion of the probe by first computing a rigid registration of all the US images together. Then we run the deformation tracking algorithm. The registration of each image of the sequence takes between 10 and 15 minutes on a standard PC running linux.

Results are presented in figure 3: on the first line, we show the original US images after the rigid registration. The second line represents the first US image deformed to match the above US image. On the last two lines, we show the MR image registered to the original US (our "ground truth") and the virtual MR produced by the tracking algorithm. To assess the quality of the tracking, we superimposed on the virtual MR images the contours of the balloons extracted from the "original" MR images.

Even if there are very few salient landmarks (all the information in the US images is located in the thick and smooth balloons boundaries, and thus the tracking problem is loosely constrained), results are globally good all along the sequence. This shows that the SSD criterion correctly captures the information at edges and that our regularized free-form deformation field is able to interpolate reasonably well in uniform areas.

When looking at the virtual MR in more details, one can however find some places where the motion is less accurately recovered: the contact between the balloons and borders of the US images. Indeed, the parameterization of the transformation and especially its smoothing are designed to approximate the behavior of a uniform elastic like body. If this assumption can be justified for the shift of brain tissues, it is less obvious for our phantom where balloons are placed into a viscous fluid. In particular, the fluid motions between the two balloons cannot be recovered. On the borders of the US images, there is virtual MR 2 virtual MR 3 virtual MR 4 virtual MR 5

Fig. 3. Tracking deformations on a phantom. In this figure, each triplet of 2D images represents 3 orthogonal views resliced from the 3D image. Top: The first 5 images of the sequence of 10 images after a rigid registration to compensate for the motion of the probe and the "virtual" US images (US 1 deformed to match the current US image) resulting from the tracking. Bottom: The "original" MR images (rigidly registered to the corresponding US images to correct for the probe motion and the phantom motion between MR acquisitions) and the virtual MR image synthetized using the deformation field computed on the US images. To assess the quality of the tracking, we superimposed the contours of the "original" MR images.

The volume of the balloons ranges from 60 to 90 ml for the ellipsoid one and 40 to 60 ml for the more complex one.

sometimes a lack of intensity information and the deformation can only be extrapolated from the smoothing of neighboring displacements. Since we are not using a precise geometrical and physical model of the observed structures like in [START_REF] Skrinjar | Real time 3D brain shift compensation[END_REF], one cannot expect this extrapolation to be very accurate.

As a conclusion from this experiment, one can say that elastic-like deformations are qualitatively well tracked in the sequence if there are some salient intensity landmarks surrounding the area of interest.

Real (pig) brain images

This dataset was obtained by Dr. Ing. V. Paul at IBMT, Fraunhofer Institute (Germany) from a pig brain at a post-lethal status. A cyst drainage has been simulated by deflating a balloon catheter with a complete volume scan at three steps. We present in figure 4 the results of the tracking. Since we have no corresponding MR image, we present on the two last lines the deformation of a grid (a virtual MR image...), to emphasize the regularity of the estimated deformation, and the deformation of a segmentation of the balloon. The registration of each image of the sequence takes between 10 and 15 minutes on a standard PC running linux.

The correspondence between the original and the virtual (i.e. deformed US 1) images is qualitatively very good. In fact, if the edges are less salient than in the phantom images, we have globally a better distribution of intensity features over the field ov view due to the speckle in these real brain images.

One should also note on the deformed grid images that the deformation found is very smooth.

To obtain a quantitative measurement of the transformation, we segmented the first image and we deformed this segmentation according to the estimated transformation field (see bottom line of Fig. 4). We can now compare the volume of the deformed balloon with its theoretical value. In fact, since the segmentation originally overestimates the balloon volume, we have to compare the ratio between the deformed volume and the original one. These images correspond to the deformation of an image of a 3D grid (a "virtual MR" image) with strips orthogonal to each 2D resliced plane: they allow to visualize the in-plane deformation for each 2D slice. Bottom: We segmented the balloon on the first image. Then, this segmentation is deformed using the transformation found and superimposed to the corresponding original US image.

The measurements indicates that we are overestimating the volume (underestimating the deformation) by 7.5% for image 2, by 3.3% for image 3, and by 30% for image 4. However, one should note that volume measurements are very sensitive as they relate to the cube of the balloon dimension: this corresponds to an error of less than one millimeter on the balloon diameter. This could be explained by an occlusion of the lower part of the balloon probably due to an air bubble trapped inside the balloon during the experience: on US 4, almost the entire lower half of the balloon is shadowed by the air bubble.

In these conditions, one cannot expect a perfect retrieval. The estimated deformation at the occlusion being computed thanks to the regularization of the deformation field from neighboring structures, it is expected to be less than the real deformations (maximal at the balloon boundaries).

Reducing the smoothing of the transformation could allow the algorithm to find a closer fit. However, this could allow some unwanted high frequency deformations due to the noise in the US images. We believe that it is better to recover the most important deformations and miss some smaller parts than trying to match exactly the images and have the possibility to create some possibly large deformations.

Discussion and conclusion

The algorithm presented here partly fills the goals of the ROBOSCOPE project: it is able to recover an important part of the deformations along the sequence and issues a smooth deformation, despite the noisy nature of the US images.

Experiments show that this allows to simulate virtual MR images qualitatively very close to the real ones. Quantitative measurements remains to be done, but it seems that an accuracy of 1 to 2 mm is achievable in the areas where there is an elastic deformation. This is encouraging since the accuracy of the clinicians without per-operative imaging is estimated to be around 3 to 5 mm. However, some improvements of the algorithm will likely be needed to cope with non-elastic deformations in the CSF, the skull, and with the introduction of the surgical tools.

We observed that the SSD criterion is well adapted for the registration of successive images in the time-sequence and performs well on our examples for the update of the global transformation. However, it is possible that other types of sequences with intensity changes may require a more complex criterion like the LCC.

The type of transformation is a very sensitive choice for such a tracking algorithm. We made the assumption of a "uniform elastic" material. This may be adequate for the brain tissues, but probably not for the ventricles and for the tracking of the surgical tools themselves. Indeed, they will penetrate into the brain without any elastic constraint with the neighboring tissues. A specific adaptation of the algorithm around the tools will likely be necessary. Another possibility for errors is the occlusion of a part of a structure visible in the US, for instance the shadowing by the endoscope.

The computation time is still far from real time for a continuous tracking of deformations during surgery but the implementation was focused on generic components in order to test different criteria and gradient descent methods. A dedicated re-implementation of the method may gain a factor 4 to 8, leading to a clinically useful tool for brain shift estimation (one estimation every minute or 2). To be further accelerated and reach real-time video-rate for instance, the algorithm must be parallelized. This would impose stronger hardware requirements but it is rather straightforward for both the computation of the image similarity and the regularization energies.

There are different parameters to tune in the algorithm but we believe that most of them could be adjusted for specific types of US images sequences. More sequences are anyway necessary to validate the estimation of the deformation.

In conclusion, we developed a tracking algorithm adapted to time sequences of US images and not only the registration of two images. Experiments on a phantom and on a real (pig) brain sequence show that the main part of the deformation is retrieved with a smooth deformation field. The image similarity criterion being independent from the type of transformation used, it could be changed in the future to better fit the assumptions on the US images depending on the application considered. We have shown here that the SSD criterion performs reasonably well in view of real-time considerations, even if a specific parallel version has to be designed in order to meet all the time requirements. However, more experiments will be needed to choose the best adapted parameterization of deformations and to validate the accuracy of the estimation.

Fig. 1 .

 1 Fig. 1. Overview of the image analysis part of the Roboscope project.

Fig. 2 .

 2 Fig. 2. The deformations computed in the tracking algorithm.

Fig. 4 .

 4 Fig. 4. Tracking deformations on a pig brain. In this figure, each triplet of 2D images represents 3 orthogonal views resliced from the 3D image. Top: The 4 images of the pig brain with a deflating balloon simulating a cyst drainage. Middle: deformation of a grid to visualize more precisely the location of the deformations found.These images correspond to the deformation of an image of a 3D grid (a "virtual MR" image) with strips orthogonal to each 2D resliced plane: they allow to visualize the in-plane deformation for each 2D slice. Bottom: We segmented the balloon on the first image. Then, this segmentation is deformed using the transformation found and superimposed to the corresponding original US image.

Acknowledgments

This work was partially supported by the EC-funded ROBOSCOPE project HC 4018, a collaboration between The Fraunhofer Institute (Germany), Fokker Control System (Netherlands), Imperial College (UK), INRIA (France), ISM-Salzburg and Kretz Technik (Austria). The authors address special thanks to Prof. Auer and his colleagues at ISM (Austria) for the acquisition of the phantom sequence, and to Dr. Ing. V. Paul at IBMT, Fraunhofer Institute (Germany) for the acquisition of the pig brain images.