Marek Bucki
email: marek.bucki@imag.fr

Olivier Palombi

Mathieu Bailet

Yohan Payan

Doppler Ultrasound Driven Biomechanical Model of the Brain for Intraoperative Brain-Shift Compensation: A Proof of Concept in Clinical Conditions

Keywords: Biomechanical modeling, Finite Element Method, Brain-shift, Neuronavigation, Doppler ultrasound

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Image guided neurosurgery is currently facing critical obstacles for large skull openings, with intraoperative brain deformations that remain hardly predictable and must be mentally taken into account by the surgeon. Deformation of brain tissues occurs in the course of surgery because of physical phenomena (dura opening, gravity, loss of cerebrospinal fluid, actions of the neurosurgeon, etc), physiological phenomena (swelling due to osmotic drugs, anesthetics, edema, etc), and other unknown causes.

The consequences of this "brain-shift" cannot be accurately anticipated preoperatively. Studies have shown that it is difficult to establish a clear correspondence between the nature of the pathology, the patient's age, the craniotomy location and the brain deformation amplitudes and direction [START_REF] Nimsky | Quantification of, visualization of, and compensation for brain shift using intraoperative magnetic resonance imaging[END_REF] [2] [START_REF] Reinges | Course of brain shift during microsurgical resection of supratentorial cerebral lesions: limits of conventional neuronavigation[END_REF]. Reinges et al. found a correlation between deformation magnitude and tumor volume, interhemispheric line shift and craniotomy dimensions, but other authors such as Nimsky et al. or Roberts et al.,[START_REF] Roberts | Intraoperative brain shift and deformation: a quantitative analysis of cortical displacement in 28 cases[END_REF], could only correlate it with the lesion volume, despite the established relationship between patient age and brain tissue mechanical properties [START_REF] Bucholz | The correction of stereotactic inaccuracy caused by brain shift using an intraoperative ultrasound device[END_REF]. Nabavi et al. [START_REF] Nabavi | Serial intraoperative magnetic resonance imaging of brain shift[END_REF] also reported that in the case of patients having undergone radiotherapy, no brain-shift was observed, probably due to tissue stiffening.

In all studies, the deformation seems to occur mainly in the direction of gravity [START_REF] Roberts | Intraoperative brain shift and deformation: a quantitative analysis of cortical displacement in 28 cases[END_REF] [1] [START_REF] Trantakis | Investigation of time-dependency of intracranial brain shift and its relation to the extent of tumor removal using intra-operative mri[END_REF]. The tissue displacement can either be a downwards sagging, due to the tissue weight coupled with lack of support from the ventricular system after cerebrospinal fluid (CSF) drainage; or an upwards bulging [START_REF] Nabavi | Serial intraoperative magnetic resonance imaging of brain shift[END_REF] which can be explained by intracranial hypertension resulting from tumor growth, called "mass effect". Nevertheless, Nabavi et al. also reported brain bulging cases without clinical indications of hypertension and Nimsky et al. [START_REF] Nimsky | Quantification of, visualization of, and compensation for brain shift using intraoperative magnetic resonance imaging[END_REF] observed bulging during interventions on meningiomas, after the tumor resection had relieved the brain from the lesion weight.

The regions affected by the brain-shift are often wider than the craniotomy and affect subcortical structures even in the case of small superficial tumors. The cortical displacements propagate deeply and can result in a compression of the ventricles. The mechanical role of the ventricles in the brain-shift process is further illustrated by the fact that the cortical and subcortical deformation amplitudes are increased when the ventricular system is opened during lesion resection, as reported by Nimsky et al. [START_REF] Nimsky | Quantification of, visualization of, and compensation for brain shift using intraoperative magnetic resonance imaging[END_REF] The deformations studied in the above mentioned works were measured using cortical landmarks or tumor margin positions. The reported magnitudes present great variability. Reinges et al. observed cortical shifts up to 15.2 mm, while Nimsky et al. reported maximal surface shifts up to 23.8 mm and deeper tumor margin displacements ranging between 7.9 and 30.9 mm [START_REF] Nimsky | Quantification of, visualization of, and compensation for brain shift using intraoperative magnetic resonance imaging[END_REF] [6] [START_REF] Reinges | Course of brain shift during microsurgical resection of supratentorial cerebral lesions: limits of conventional neuronavigation[END_REF]. Nabavi et al. measured maximal shifts up to 5 cm in the resection areas.

With regards to the dynamic aspect of brain-shift, all authors agree that the deformation occurs quickly after dura opening [START_REF] Trantakis | Investigation of time-dependency of intracranial brain shift and its relation to the extent of tumor removal using intra-operative mri[END_REF] [START_REF] Reinges | Course of brain shift during microsurgical resection of supratentorial cerebral lesions: limits of conventional neuronavigation[END_REF] and evolves throughout the intervention, which makes it necessary to monitor the brain deformation as often as possible in order to accurately capture the cortical and subcortical changes [START_REF] Nabavi | Serial intraoperative magnetic resonance imaging of brain shift[END_REF].

At the end of the intervention, after tumor removal, Nabavi et al. observed brain bulging and volume increase in all studied cases, which the authors attribute to edema or expansion of previously compressed tissues.

To our knowledge, the complex brain-shift issue is not addressed by any of the commercially available neuronavigators and the problem remains unresolved. As a consequence, routinely used image guided neuronavigators can be safely used only at the preliminary stage to precisely define the craniotomy location, or for interventions where only a small borehole needs to be drilled in the patient's skull, which strongly limits CSF leakage and hence brain-shift [START_REF] Nimsky | Quantification of, visualization of, and compensation for brain shift using intraoperative magnetic resonance imaging[END_REF]. In most cases, however, the divergence between pre-operative data and actual patient situation puts the surgical procedure accuracy at stake.

In order to compensate for the brain-shift, medical images of the patient's head can be reacquired in the operating room (OR) using intraoperative Magnetic Resonance Imaging (iMRI) [START_REF] Ferrant | Registration of 3-d intraoperative mr images of the brain using a finite-element biomechanical model[END_REF] [8] [START_REF] Chu | Intraoperative magnetic resonance imaging-guided neurosurgery[END_REF] [10] [11] [12]. This approach suffers from many limitations. First, due to its high cost, only a few clinical centers have access to this technology. Second, from the practical point of view, long patient installation and image acquisition time of a few minutes make it difficult to acquire a sufficient number of images needed to accurately track the evolution of the brain-shift during the intervention [START_REF] Wirtz | Image-guided neurosurgery with intraoperative mri: update of frameless sterotaxy and radicality control[END_REF] [START_REF] Steinmeier | Intraoperative magnetic resonance imaging with the magnetom open scanner: concepts, neurosurgical indications and procedures -a preliminary report[END_REF]. Finally, the resolution and scope of the produced images are usually smaller than those of the available pre-operative data and, as a consequence, information is lost.

A number of solutions have been proposed in the literature in order to intraoperatively update, and take advantage of the planning and other pre-operatively available data. The earliest proposed algorithms deformed the pre-operatively acquired images using image-based models. Different non-rigid registration methods were therefore provided to match intraoperative images (mainly MRI exams) with pre-operative ones [START_REF] Hata | Three-dimensional optical flow method for measurement of volumetric brain deformation from intraoperative MR images[END_REF] [16] [START_REF] Shattuck | BrainSuite: an automated cortical surface identification tool[END_REF]. More recently, biomechanical models of the brain tissues were proposed for constraining the image registration by inferring a volumetric deformation field from the correspondences between contours [START_REF] Kyriacou | Nonlinear elastic registration of brain images with tumor pathology using a biomechanical model[END_REF] [19] and/or surfaces [START_REF] Ferrant | Serial registration of intraoperative MR images of the brain[END_REF] in the images to be registered.

In their arguments against the exorbitant cost of intraoperative MRI imaging devices, some authors have proposed to couple the biomechanical model of the brain with low-cost readily available intraoperative data [START_REF] Dumpuri | An atlas-based method to compensate for brain shift: Preliminary results[END_REF] such as cortical surfaces recovered from laser-range scanner systems [START_REF] Audette | An integrated range-sensing, segmentation and registration framework for the characterization of intra-surgical brain deformations in image-guided surgery[END_REF] [START_REF] Miga | Cortical surface registration for image-guided neurosurgery using laser-range scanning[END_REF] or intraoperative ultrasound [START_REF] Comeau | Intraoperative ultrasound for guidance and tissue shift correction in image-guided neurosurgery[END_REF]. Given the high cost of iMRI devices, this proposal seems appealing from a very practical point of view; however, it gives the biomechanical model a crucial and very central position. This solution requires that a strong modeling effort is carried out during the design of the brain model and its validation against clinical data.

An interesting approach was proposed and evaluated by Reinertsen et al. in [START_REF] Reinertsen | Clinical validation of vesselbased registration for correction of brain-shift[END_REF], where the authors suggest the use of intraoperative 3D Doppler US imaging to track the brain vascular tree and apply a physical Thin Plate Splines (TPS) deformation to the pre-operative data in order to compensate for the brain-shift. Doppler US imaging of the brain during neurosurgery seems to be a promising approach as tumor growth often leads to vascular remodeling and creation of neo-vessels derived from the vascular tree to bring nutriments to the cancer cells [START_REF] Zagzag | Neovascularization and tumor growth in the rabbit brain[END_REF] [START_REF] Hjelmeland | Twisted tango: brain tumor neurovascular interactions[END_REF]. In those cases the new vessels at the vicinity of the tumor thus constitute excellent natural fiducials which can be localized to track the tumor margin. Due to its ease of use, US imaging can be repeatedly resorted to for performing updates of the pre-operative data whenever the brain structure has undergone major changes.

In this article, we propose to compute the deformation of a patient specific Finite Element (FE) biomechanical model of the brain using the sparse displacements field intraoperatively recovered from Doppler US vascular tree tracking. The FE modeling framework allows us to more easily take into account complex boundary conditions such as the contact between the cortex and the skull. Furthermore, the physical properties of the tissues such as the stiffness difference between gray and white matter and the ventricles can be modeled. Through the imaging of the patient's brain by Diffusion Tensor Imaging it is also possible to take into account the strong anisotropy of the white fiber tracts within the parenchyma [START_REF] Clatz | Realistic simulation of the 3-d growth of brain tumors in mr images coupling diffusion with biomechanical deformation[END_REF]. Unlike the bending energy implemented in TPS the physical parameters of a FE biomechanical model are easier to relate to the actual mechanical properties of the modeled tissues and in our opinion, this approach is better suited for addressing the problem of brain-shift.

As a low-cost solution was sought and 3D Doppler ultrasound was actually not available in our operating theater, the study was carried out using standard 2D Doppler imaging. This restriction raised issues related to freehand US data acquisition, such as noise and presence of outliers in the processed data, as discussed below. The technique described in this article was developed with a special focus on system reactivity and robustness. A clinical evaluation showed the ability of the neuronavigator for dealing with the demanding conditions of use in the OR.

This article is organized as follows. Section 2 gives an overview of the system, and its technical aspects are further detailed in §3 for the pre-operative phase, and in §4 for the intraoperative phase. Section 5 illustrates the neuronavigator capabilities on a clinical case and finally, in §6 the results are discussed and some perspectives of future works are outlined.

System overview

The system workflow is summarized in Fig. 1. It can be divided into a preoperative and an intraoperative phase. The pre-operative phase is represented by light gray boxes in the figure and the intraoperative phase is shown as dark gray boxes.

Before the intervention, MR images of the patient's brain are taken and the resulting volume is used to define the surgical strategy, or "planning". The neuronavigation is performed using an image guided binocular microscope, the Surgiscope system4 . During surgery the position of the microscope is tracked with respect to the preoperative patient MR images. Due to the brain-shift effect however this neuronavigator is only used to position the craniotomy and once an important shift has been observed the navigation feature is abandoned by the surgeon.

In addition to routine neurosurgical planning, the pre-operative data are processed by our system in order to:

-Recover the brain vascular tree initial configuration from Magnetic Resonance Angiography (MRA) data ( §3.1); -Build a patient-specific biomechanical FE model for brain-shift simulation ( §3.2).

The MRA sequence used to segment the preoperative vascular tree is not used routinely for neuronavigation and was acquired specifically for this study.

Once the patient is installed in the OR with the head firmly fixed in a steady frame ( §4.1), a rigid registration procedure is carried out to recover the position of the pre-operative images with respect to the physical position of the head in the OR. The coordinates in image space of the segmented brain vascular tree are converted into a patient related metric frame, called "head reference" ( §4.2). The resulting initial configuration of the vascular tree is referred to as "reference vascular tree".

During the intervention, whenever the effects of the brain-shift make it necessary, the surgeon performs a localized 2D Doppler US sweep of the region of interest. The spatially scattered 2D US images form a so-called "pseudo-volume" which segmentation makes it possible to capture a 3D shape of the deformed vascular tree ( §4.3).

Once the deformed vascular tree has been elastically registered to the reference vascular tree using a robust non-rigid registration procedure [START_REF] Bucki | A fast and robust patient specific finite element mesh registration technique: application to 60 clinical cases[END_REF], the out-put displacements field is processed to remove outliers resulting from erroneous segmentation or registration ( §4.4). The vectors describing the vessel displacements are defined to stem from the reference vessel centerlines and point to the corresponding centerlines in the deformed vascular tree configuration.

Next, the biomechanical model is used to smooth the sparse displacements field and compute the overall deformation of the brain. The FE deformation computation is optimized in order to provide to the surgeon a relevant representation of the organ within seconds after the US images acquisition ( §4.5).

Finally, the pre-operative data are updated according to the FE model deformation and displayed to the surgeon who performs a verification of the estimated brain-shift using Doppler US imaging ( §4.6). If a satisfactory accuracy is reached in the region of interest, the planning update proposed by the system is accepted and the intervention resumes relying on a more faithful representation of the patient's brain.

Once the planning has been accepted by the surgeon the therapeutical act is performed, which usually involves tissue resection. In order to accurately assist the clinician throughout the intervention a neuronavigation system must be able to take into account the tumor removal and compensate for the additional shift that the resection induces. A novel approach to intraoperative modeling of tissue resection is proposed and some preliminary results are presented ( §4.7).

Pre-operative phase

In this section we describe the data processing steps performed before the intervention: identification of the reference vascular tree from 3D Time-of-Flight Magnetic Resonance Angiography (3DToF MRA), and construction of a patient specific FE biomechanical model.

Magnetic resonance angiography acquisition and analysis

In the context of brain tumor removal surgery, Magnetic Resonance Imaging is the preferred modality chosen for the high level of detail among soft tissues that it is able to distinguish. T1 and T2 image processing sequences strongly contrast tumors or other lesions with normal tissues and make it possible to define an accurate neurosurgical planning. On a T2-weighted scan, tissues containing fat, water and fluid appear as bright hypersignal. As damaged tissues tend to develop edema, T2-weighted sequences are sensitive to pathology and are usually able to distinguish pathologic from normal tissues.

Vascular tree imaging using 3DToF MRA represents one of the most interesting diagnostic MR capabilities in the neurovascular field. In order to better visualize the distal portions of the brain vascular tree that usually exhibit slower flow rate, paramagnetic gadolinium-based contrast agents can be used [START_REF] Yano | Gadolinium-enhanced 3D timeof-flight MR angiography. Experimental and clinical evaluation[END_REF]. This enhancement is relevant here as tumor margin shifts are mostly revealed by the position of small neo-vessels. Unlike Computed Tomography (CT), 3DToF MRA is a non invasive 3D angiography technique as there is no ionizing radiation involved and, to our knowledge, no significant side effects of radio waves or magnetic fields on the human body have been reported to date.

3DToF MRA is therefore a well suited pre-operative modality as it is riskfree and able to reveal important arteries as well as finer vessels near the brain lesion. For patients who do not require an MRA exam before the intervention the tradeoff associated with the extra imaging cost is a more accurate surgical procedure performed under the supervision of a system capable of compensating for a possible brain-shift.

Given the high contrast quality of the images produced by 3DToF MRA, a basic semi-automatic segmentation technique was used to identify the vascular tree, relying on simple gray level thresholding [START_REF] Reinertsen | Vessel driven correction of brain shift[END_REF]. Ultimately this procedure should be fully automatic as the clinicians can hardly dedicate time to perform this kind of pre-processing task. In our study this work was not done by the surgeon but by an expert assistant.

The segmentation procedure was performed in four steps, as follows:

1. A gray level threshold value was manually adjusted so as to visualize a maximal number of vessels while limiting the inclusion of other, non-vascular structures; 2. As the diameter of larger vessels changes throughout the heart beat cycle due to the arterial wall compliance, only the vessel centerlines were retained as fiducials whose position is common in both the MR and US imaging modalities. The centerlines were defined as mass centers of the thresholded areas; 3. Artifacts were filtered using a proximity criterion. From the set of segmented centerlines those further away than 10 mm from other centerlines were removed which discarded most of the isolated segmented features with a lesser probability to pertain to the vascular tree; 4. Finally, the operator checked the consistency of the vascular tree with the possibility to manually add or remove vessel centerlines.

Fig. 5-c shows the points cloud formed by the segmented vessel centerlines obtained during the evaluation of our system ( §5). The resolution of the 512 × 512×120 MRA volume was 0.39×0.39×0.7 mm. The resulting 2600 points set was dense enough so that continuity of the vascular structure could be visualized. The exact connectivity however needed not to be computed as the elastic registration procedure used in this study is able to deal with both connected and sparse structures. Besides, given the high density of information retrieved in the region of interest this process was considered unnecessary.

Patient specific FE model

Biomechanical modeling of the brain raises two issues. First, an appropriate model of the biomechanical behavior of the tissues under scrutiny must be implemented using partial derivative equations (PDEs) formulated in the framework of Continuum Mechanics and second, a faithful geometrical representation of the modeled organ must be achieved in order to compute an accurate response of the model under the prescribed load. The deformation field, solved with PDEs is usually estimated using the well known FE method [START_REF] Belytschko | Nonlinear Finite Elements for Continua and Structures[END_REF] which relies on a discretization, or "mesh", of the modeled domain.

The biomechanics of the brain soft tissues have been widely studied both in-vitro and in-vivo, and their non-linear and nearly incompressible nature has been established [START_REF] Metz | A comparison of the elasticity of live, dead, and fixed brain tissue[END_REF] [33] [START_REF] Miller | Biomechanics of brain for computer integrated surgery[END_REF] [35] [START_REF] Soza | Estimating mechanical brain tissue properties with simulation and registration[END_REF]. A number of potential elastic energy formulations have been proposed by the authors nevertheless, despite recent improvements in that area [START_REF] Schiavone | In vivo measurement of human brain elasticity using a light aspiration device[END_REF], the rheological parameters remain difficult to assess in-vivo and intraoperatively. As a result of that uncertainty a great variety of values for the Young's modulus can be found in the literature: E = 8196P a in [START_REF] Soza | Estimating mechanical brain tissue properties with simulation and registration[END_REF], E = 35kP a in [START_REF] Metz | A comparison of the elasticity of live, dead, and fixed brain tissue[END_REF] or E = 313kP a in [START_REF] Eckabert | Brain shift correction based on a boundary element biomechanical model with different material properties[END_REF]. The recent in-vivo studies by Miller et al. and Schiavone et al. seem to find an agreement for a value of E around 1kP a.

As system robustness and short computational times are crucial during surgery, a linear model of the brain tissues was chosen here and the hypothesis of both small deformations (mechanical linearity) and small displacements (geometrical linearity) was adopted. Clatz et al. [START_REF] Clatz | Brain tumor growth simulation[END_REF] have shown that under the assumption that the brain deformation does not exceed 10% a linear model of the parenchyma with E par = 694P a and ν par = 0.4 does not diverge from the nonlinear law proposed in [START_REF] Miller | Biomechanics of brain for computer integrated surgery[END_REF] by more than 7%. The high anisotropy of white matter was ignored. Given its free circulation within the ventricles, the cerebrospinal fluid (CSF) was modeled as a very compressible medium, taking the values E csf = 10P a and ν csf = 0.05 as suggested by Clatz et al. It is important to note that since the biomechanical model is displacements-driven and implements a linear material only the relative value of E par and E csf is relevant. One of the moduli could be scaled to 1 without affecting the produced deformations.

The biomechanical parameters used in this study, E par = 694P a, ν par = 0.4, E csf = 10P a and ν csf = 0.05, were taken from prior works in the literature [START_REF] Clatz | Brain tumor growth simulation[END_REF] and are not patient specific. On the other hand, the geometry of the brain, the shape and position of the ventricles as well as the boundary conditions resulting from the planned position of the craniotomy were adapted to each specific clinical case, as described below.

The brain parenchyma and hypointense ventricles were both segmented using the 3D Slicer5 software from the T1-MR volume acquired for the routine neuronavigation procedure. A simple region growing algorithm was applied after a seed was manually placed in the corresponding brain area. The automatically segmented binary volume was cleaned up using morphological opening and closure operators and finally the largest connected component was retained for each considered structure.

The FE mesh was generated using the MMRep FE mesh registration technique [START_REF] Bucki | A fast and robust patient specific finite element mesh registration technique: application to 60 clinical cases[END_REF] that deforms a template, or "Atlas", brain model onto the cortical surface of the patient's brain segmented from the preoperative MR dataset. The brain shape representation accuracy was measured as the mean distance between the mesh surface nodes and the corresponding segmented surface. A submillimetric accuracy was achieved producing a 2804 nodes and 14322 elements tetrahedral patient specific mesh. As the ventricular system has not yet been defined in the Atlas mesh used here, instead of being registered simultaneously with the cortical surface the segmented ventricles were modeled by implementing the values for E csf and ν csf for the mesh elements fully included in their volume.

The boundary conditions were imposed on the patient specific model according to the planning defined by the surgeon. All interior mesh nodes had unconstrained displacements and boundary conditions were only defined for cortex mesh nodes.

As the orientation of the patient's head in the operating room with respect to the direction of gravity was known beforehand, the brain-skull contact boundary conditions were set as follows. The contact between the skull and the mesh nodes on the lower half of the cortical surface was modeled as a frictionless sliding on the tangent plane at the considered skull point. The nodes on the brain stem were fixed as suggested in [START_REF] Wittek | Patient-specific model of brain deformation: application to medical image registration[END_REF] and the influence of gravity was ignored as its influence was to be revealed by the observed displacement of the vascular tree. Finally all other mesh nodes were set to move freely in all directions.

The deformation precomputations were carried out before the intervention on the produced patient specific biomechanical model so as to take advantage of the linear formulation [START_REF] Cotin | Real-time elastic deformations of soft tissues for surgery simulation[END_REF], as described in further detail in §4.5.

Intraoperative phase

System set-up in the operating room

The system layout in the operating room (OR) is illustrated in Fig. 2. A Polaris6 passive optical localizer is used to track the position of the surgical tools and the US transducer in the OR space during surgery. The localization error of the device has an RMS of 0.35 mm. As the localizer needs a clear line of sight to the operating field, it is placed between the surgeon and the assistant, and pointed towards the patient's head.

The "head" referential is defined by a rigid body attached to the head clamp holding the patient's head (Fig. 2). During the US scanning of the brain, the transducer is localized by means of a rigid body attached to it, defining the "US" referential. The transducer's rigid body and ultrasound cable are isolated from the operating field using a sterile endoscope cover. The rigid body is designed so that the sterile reflective markers7 can be clipped on from the outside through the plastic sheet to avoid occlusions or unwanted reflections from the sterile cover.

The localizer is connected to the neuronavigator ("PC" in the figure) which continuously retrieves the relative position of all tracked devices. The ultrasound

Vascular tree rigid registration

The first step taken at the beginning of the intervention is the conversion of the pre-operatively segmented vascular tree coordinates from image space into the patient's physical space, measured with respect to a "head" referential rigidly attached to the patient's head through a head clamp and materialized by the "head" rigid body. All subsequent position computations are performed in the "head" referential.

Prior to the Surgiscope MR image acquisition, a set of 7 to 10 skin markers are placed on the patient's head. These markers are clearly visible in the MRI and their image coordinates are manually identified by the surgeon during the pre-operative planning phase. Once the patient's head is fixed in the head clamp and the "head" rigid body is mounted upon it, the landmark points are localized in the patient's space by placing the tip of a calibrated pointer at the center of each fiducial and recording the position with respect to the "head" referential. The resulting fiducial points in patient's physical space are rigidly registered on their counterparts in MRI image space using an Arun least squares minimization [START_REF] Arun | Least-squares fitting of two 3-d point sets[END_REF].

Doppler US acquisition and analysis

The Doppler analysis of the ultrasound signal reveals functional information about the organ as it measures fluid flow velocity and direction. Color Doppler ultrasound is mostly used to measure the rate of blood flow through the heart and major arteries. The color coded flow areas are overlaid on the anatomical brilliance mode, or "B-mode" images to form a so-called "M-mode" (multi-mode) composite image. Power Doppler is an enhancement of the technique that is up to 5 times more sensitive in detecting blood flow than standard color Doppler. In power Doppler mode smaller distal vessels can thus be visualized more accurately. Unlike in Color Doppler where the blood flow direction are encoded in red or blue depending on the orientation of the flow with respect to the US slice, in power Doppler mode the flow orientation cannot be visualized and the flow speed alone is represented usually in shades of orange.

Power Doppler tracking of the brain vessels in the context of neurosurgery has already shown promising results. A vessel driven approach has been developed by [START_REF] Lee | Vascular tree extraction from mra and power doppler us image volumes[END_REF] and more recently a clinical validation of a similar procedure has been presented in [START_REF] Reinertsen | Validation of vessel-based registration for correction of brain shift[END_REF]. The study carried out on 5 patients by Reinertsen et al. showed that a TPS deformation coupled with intraoperative Doppler imaging was able to correct the brain shift deformation to an accuracy within 1.25 mm in about 60 seconds. Furthermore, a recent study on the 3D reconstruction of the vascular tree using intra-cranial power Doppler ultrasound concluded that 'the anatomical relation between the vascular structures and the tumor was clear ' [START_REF] Yasuda | Intraoperative three-dimensional reconstruction of power doppler vascular images[END_REF], which shows that the displacement information yielded by brain vessels tracking can be used to effectively recover the position of a tumor during surgery. In the remainder of the article we will simply refer to power Doppler imaging as "Doppler".

Before the intervention, a rigid body is firmly attached to the US transducer and the ultrasound probe is calibrated using the freehand membrane scan technique [START_REF] Lango | Ultrasound guided surgery: Image processing and navigation[END_REF]. This time-consuming procedure is carried out beforehand and the resulting matrix is saved. At the beginning of the intervention the calibration information is quickly recovered at neuronavigator start-up.

Whenever a large deformation occurs and the pre-operative data no longer reflect the actual organ configuration, the surgeon can perform an ultrasound sweep of the region of interest. To this end the craniotomy cavity is filled with saline solution routinely used to clean the surgical wound throughout the surgery. This fluid is not only an excellent coupling medium but it also prevents the probe from touching the cortical surface, which could create unwanted tissue deformation and induce bias in the brain-shift estimation. Of course this technique assumes that the craniotomy is facing up and can be filled with fluid.

To acquire the set of Doppler US images, the hand-held probe is rotated by the surgeon so as to sweep the largest possible volume around the region of interest while the US transducer is being localized with respect to the patient's head reference. The set of recorded images along with their spatial positions define a so-called "2.5D" or "pseudo-volume".

After the set of images has been acquired, the Doppler US pseudo-volume is segmented as follows. First, the colored pixels are isolated from the B-mode background by a thresholding performed on the saturation value computed after converting the Red-Green-Blue (RGB) signal into the Hue-Saturation-Value (HSV) color space. The retained pixels are those with a saturation greater than 5%, a manually adjusted value most likely hardware-dependent.

Once the Doppler signal is extracted from the images, a shape analysis is performed on the sets of connected pixels, or "stains", corresponding to the brain vessels. First, all images containing stains with an area greater than 10 mm 2 are discarded from the segmentation. This step eliminates most of the Doppler artifacts due to excessively fast freehand motion that appear in the images as large colored stains spreading radially, along the US waves propagation direction. In a second selection step only the stains having an aspect ratio between 1 and 2 are retained, which corresponds to a ±30 • maximal angle between the US image plane and an ideal cylindrical vessel axis. Finally, the centers of the selected stains are converted into the patient's space physical coordinates using probe calibration and position information. The final result is a points cloud describing the deformed vessel centerlines, as shown in Fig. 6.

The segmentation procedure described above discards most of the noise inherent to the Doppler modality, yet signal due to motion of fluids other than blood can also be recorded. In our experience cerebrospinal fluid (CSF) circulation within the ventricles could be observed and produced round, vessel-like Doppler stains. In those cases the ventricular origin of the Doppler signal was established after registering the recorded US images with pre-operative MR scans of the patients.

Hardly any further information on the intraoperative vascular tree structure can be recovered from a 2D Doppler US acquisition. This is mainly due to the fact that the patient's heart palpitation generates a varying blood flow in the brain arteries causing small vessels to appear in one slice and be absent from a slice immediately adjacent in space and time, depending on the ultrasound sensitivity and probe frequency. The large deeper arteries, such as those in the Circle of Willis, are present in all the Doppler US slices independent of the heartbeat but the sparsity of the data describing the thinner distal vessels does not allow a proximity search among the gathered centerlines to establish connectivity relationships in the outer cortical areas. The nature of the intraoperative data set describing a portion of the deformed brain vascular tree is thus a sparse cloud of points. Although structural information is missing, the tracking of the vascular tree yields a great quantity of valuable information widely spread throughout the region of interest.

Vascular tree elastic registration

The data gathered from the pre-and intraoperative brain vascular tree identification steps described above are quite different. The pre-operative vascular tree shape and connectivity is almost entirely recovered with little noise on the signal, the MRA modality encompasses almost the entire volume of the patient's head, and the segmentation procedure discards almost all outliers. The typical voxel size is approximately 0.4×0.4×1.0 mm. The Doppler US images recorded during the manual sweep, on the other hand, only cover a portion of the vascular tree volume. Depending on the op-erator's skills the pseudo-volume slices can be more or less evenly spaced with possible multiple sweeps over a given region while leaving out other areas. The error on the position of the intraoperatively computed brain vessels is the sum of the errors on the spatial probe localization and the uncertainty on the probe calibration itself. Furthermore, a high number of outliers can be found in the data mostly due to soft tissue displacements likely to occur during the sweep, unwanted contact of the hand-held probe with the cortex or detection of fluid flows other than blood in the brain arteries. The pixel resolution is finer than that of the MRA, with a typical US image pixel size about 0.2×0.2 mm.

The task of the elastic registration algorithm is thus to establish a correspondence between the sparse, incomplete, redundant and noisy intraoperative "current" configuration of the vascular tree and the clean, complete pre-operative "reference" representation of the same structure. From this correspondence the set of displacement vectors stemming from the reference will be derived which point to the current, intraoperative configuration.

Our team has recently proposed a fast and robust registration technique estimating a correspondence between the so-called "source" and "destination" data sets, noted S and D respectively. The resulting diffeomorphic deformation R is iteratively assembled on a hierarchical regular grid and minimizes a registration energy E(R) defined as the sum over the source points in S of the Euclidean distances to the closest point in the destination set D, i.e.:

E(R) = s∈S d(R(s), D) (1)
At each elastic registration iteration, space distortion is limited by locally monitoring the potential elastic energy of a virtual solid enclosing the processed data and choosing the deformation leading to the smallest elastic potential increase. This technique was successfully used in the context of automatic FE mesh generation to adapt generic, or "Atlas" FE meshes to patient-specific organ shapes recovered from pre-operative or intraoperative imaging [START_REF] Bucki | A fast and robust patient specific finite element mesh registration technique: application to 60 clinical cases[END_REF].

The asymmetrical definition of E(R) makes it possible to deal with situations where S defines a subset of the structure featured in the destination set D. In our case, as the intraoperatively recovered vascular tree is a portion of the pre-operatively identified structure, the source points set S is obtained by the segmentation of the Doppler US images, while D is the set of vessel centerlines yielded by the segmentation of the pre-operative MRA.

Thanks to the mechanical regularization, the elastic registration procedure is resistant to noise; nevertheless, the presence of outliers in S, as discussed above, may lead to associations between source and destination points that do not correspond to reality. It is therefore important to remove from the produced displacements field these mis-registrations before proceeding to the biomechanical model deformation. To this end, the post-processing filter described below relies on the computation of the inverse registration function, which can be achieved in our framework with an arbitrary accuracy. Source and destination point associations that do not meet a given accuracy standard are removed in two steps, by applying a "forward" filtering of the displacements field followed by a "backwards" filtering. Let d max be the maximum acceptable distance between a registered source point and its closest counterpart within the destination points set D. This value was set to 1 mm in our implementation. The forward filtering step eliminates the points in S that did not reach the destination set D within accuracy d max . The forward-filtered set of source points is thus:

S f = {s ∈ S | d(R(s), D) < d max } (2)
The goal of the registration procedure, as stated above, is to generate a displacements field reflecting the movement of the reference vascular tree, defined by D, towards its intraoperatively observed shape S. The backwards filter identifies those points in D that reach by inverse registration a vessel centerline in S f within accuracy d max . To this end, the inverse R -1 of the registration function is computed and the backwards-filtered set of destination points is defined as:

D b = {d ∈ D | d(R -1 (d), S f ) < d max } (3)
Finally, the input to the biomechanical model deformation procedure described in the next section is the set of vessel displacement V, built upon D b , and defined as:

V = { - → v = -→ ab | a ∈ D b , b = R -1 (a)} (4)

Biomechanical model deformation

The displacements field given by the elastic registration procedure has the following characteristics:

-The origin of each displacement vector is not necessarily located at a FE mesh node as the position of the reference vascular tree was not taken into account during the patient specific FE mesh generation; -The estimation of the vascular tree displacement by the elastic registration procedure relies on a Euclidean similarity measure between the pre-and intraoperative point clouds and therefore ignores the biomechanical properties of the underlying tissues.

The first point makes it difficult to rely on classical boundary conditions imposition techniques such as the penalty method, the Lagrange-Multipliers [START_REF] Hughes | The finite element method: Linear static and dynamic finite element analysis[END_REF], the mesh condensation or sub-structuration proposed in [START_REF] Kardestuncer | Finite element handbook[END_REF] [START_REF] Bro-Nielsen | Finite element modeling in surgery simulation[END_REF] where the model is deformed by applying a set of displacements at the mesh nodes. The second aspect raises the issue of the integration of uncertain vessel displacements within a biomechanical model which should not only extend the partial displacement information to the entire organ volume, but also provide a way of smoothing the displacements vector field. The following 3-steps algorithm has been implemented to address these constraints.

1. Let E d be the set of indices of mesh elements containing at least one displaced vessel point. For each element e in E d the displacements of its interior points are converted into nodal displacements and stored as a nodal displacement vector D e ; 2. For each element e in E d the impact of D e on all the mesh nodes is computed and stored as a global nodal displacement vector U e ; 3. The brain model deformation U is finally computed as the linear combination of all U e , e ∈ E d , that minimizes the squared errors on the vessel displacements V defined in Eq. 4.

The three steps are detailed below.

1 -Conversion of inner displacements into nodal displacements D e . In order to convert the displacement of interior points within a mesh element into nodal displacements, the technique proposed in [START_REF] Hsu | Direct manipulation of free-form deformations[END_REF] [50] was used. This approach leads to a solution with a least squares error on the imposed inner displacements when the problem is over-determined, and a set of nodal displacements with a minimal magnitude when the problem is under-determined. Let's consider an element e in E d , and let

V e ⊂ V, V e = { - → v i = --→ a i b i } i=1,.
..,Ie be the set of I e ≥ 1 displacements applied to the interior points {a i } i=1,...,Ie in element e. The value of each shape function ϕ n associated to the n th node of the tetrahedron e can be computed for each vessel point a i as ϕ n i = ϕ n (a i ), n = 1, . . . , 4, i = 1, . . . , I e .

The searched nodal displacement vector D e has 12 rows and can be written as D e = (d

→ v i = (v i1 , v i2 , v i3
) T for all i = 1, . . . , I e , then D e is the solution of the linear system:

           ϕ 1 1 0 0 ϕ 4 1 0 0 0 ϕ 1 1 0 • • • 0 ϕ 4 1 0 0 0 ϕ 1 1 0 0 ϕ 4 1 . . . . . . . . . ϕ 1 Ie 0 0 ϕ 4 Ie 0 0 0 ϕ 1 Ie 0 • • • 0 ϕ 4 Ie 0 0 0 ϕ 1 Ie 0 0 ϕ 4 Ie                       d 1 1 d 1 2 d 1 3 . . . d 4 1 d 4 2 d 4 3            =            v 11 v 12 v 13 . . . v Ie1 v Ie2 v Ie3            (5)
Or, in a shorter form: H e D e = V e . The dimensions of the matrix H e are 3I e × 12 and thus the system is under-determined if 1 < I e ≤ 4 and over-determined if I e > 4. In both cases the unknown vector D e is found using the pseudoinverse of the matrix H e , noted H + e and defined in the under-determined case as:

H + e = (H T e H e ) -1
H T e , and in the over-determined case as: H + e = H T e (H e H T e ) -1 . The searched nodal displacement vector is thus:

D e = H + e V e (6)
2 -Computation of the overall model deformation U e . The fast computation of the overall model response to a set of imposed nodal displacements D e estimated for element e at the previous step relies on the principle of superposition, or linear combination of pre-computed deformations [START_REF] Bro-Nielsen | Finite element modeling in surgery simulation[END_REF].

In order to compute the displacements U e of all the mesh nodes resulting from the displacement D e of the nodes in tetrahedron e, the mesh nodes are re-numbered so that U e can be written as a compound vector formed by the known nodal displacements vector D e and the unknown displacements vector U e , thus:

U e = D e U e (7)
In the absence of external nodal forces the linear system defining the mechanical behavior of the model can be written using the reordered stiffness matrix K e as K e U e = 0. The matrix K e can be further decomposed, leading to:

R e S T e S e T e D e U e = 0 0 (8)
The lower part of the above linear equation gives the expression of the unknown nodal displacements as function of the imposed displacements:

U e = -T -1
e S e D e . Once U e computed and the overall mesh deformation U e assembled, its rows are permuted so as to restore the mesh nodes initial order.

The inversion of each matrix T e , e ∈ E d , is computationally intensive and in order to comply with the intraoperative time constraints each solution U e is assembled as a linear combination of precomputed "elementary" deformations. To this end, let's consider the decomposition of the vector D e on a canonic base:

D e = 3 i=1 4 n=1 d n i δ n i ( 9
)
where the coordinates of δ n i are all 0 except the (3n + i) th which is equal to 1. The expression of the searched nodal displacements vector U e thus becomes:

U e = 3 i=1 4 n=1 d n i -T -1 e S e δ n i = 3 i=1 4 n=1 d n i U n ei (10)
Each solution U e can be intraoperatively assembled within a fraction of second assuming that the set of vectors {U n ei } e,i,n has been precomputed for all the elementary displacements δ n i , in all dimensions i = 1, . . . , 3, for all nodes n = 1, . . . , 4, and within all elements e ∈ E d involved in the deformation of the patient's vascular tree.

These precomputations are patient-specific as they depend on the shape of the brain FE mesh and the position of the reference vascular tree identified from the pre-operative MR images. Although computational performance is not a hard constraint at the pre-operative stage, the time required for the inversions of the many {T e } e matrices was significantly reduced by resorting to the Conjugate Gradient iterative solution scheme [START_REF] Shewchuk | An introduction to the conjugate gradient method without the agonizing pain[END_REF]. Performance was further enhanced using a simple Jacobi preconditioning. More details on the precomputation of the elementary deformations are given in §5.1.

-Assembly of the overall brain model deformation U.

As many elements among E d share nodes, their mutual influence must be taken into account when assembling the global brain model deformation U from the set of individual deformations {U e } e∈E d computed at the previous step. An exact solution to the superposition of a set of nodal displacements has been proposed in [START_REF] Cotin | Real-time elastic deformations of soft tissues for surgery simulation[END_REF] and requires the solution of a 3m × 3m linear system, where m is the number of displacements imposed on the FE mesh nodes. This technique, proposed in the surgical simulation context, is well-suited when m is small, for example, when modeling the interaction of a surgical tool with the surface of an organ through a limited number of mesh nodes. In our case, however, the number m of nodes with imposed displacements is usually large (see §5.2) and for better intraoperative system reactivity this computational complexity must be reduced. The proposed technique is described below.

The assembly of the final brain model deformation U from the elementary mesh deformations {U e } e∈E d is driven by the minimization of the squared error on the imposed vessel displacement field V (Eq. 4). The number of degrees of freedom to be determined is reduced as only one single scalar weight ρ e affects each considered deformation vector U e . Let N E be the number of elements in E d and ρ = (ρ 1 , . . . , ρ N E ) be the vector of scalar weights. The expression of the nodal displacements vector U becomes:

U = U(ρ) = e∈E d ρ e U e (11)
If we index the set of vessel displacements V defined in Eq. 4 so that V = { -→ v i } i∈I , then the following squared error term E(ρ) must be minimized:

E(ρ) = i∈I - → w i (ρ) -- → v i 2 (12)
where -→ v i = --→ a i b i is the vessel displacement estimated by the elastic registration and -→ w i (ρ) is the displacement produced at point a i by the FE model deformation (Eq. 11). Each vector -→ w i (ρ) can be rewritten as:

- → w i (ρ) = n∈Ni ϕ n i   U[x(n)] U[y(n)] U[z(n)]   = e∈E d ρ e n∈Ni ϕ n i   U e [x(n)] U e [y(n)] U e [z(n)]   = e∈E d ρ e -→ W ei (13
) where N i is the set of mesh node indices affecting the position of point a i , ϕ n i = ϕ (a i ) is the value of the weight function associated to node n and evaluated at a i , and x(n), y(n) and z(n) are the indices of the displacement of node n within an overall nodal displacement vector U or U e (the square braces denote the indexation operator).

Using the above identity, the expression given in Eq. 12 can be rewritten as:

E(ρ) = i∈I e∈E d ρ e -→ W ei • e∈E d ρ e -→ W ei -2 - → v i • e∈E d ρ e -→ W ei + - → v i 2 ( 14
) and can further be transformed into a quadratic form:

E(ρ) = ρ T Aρ -2ρ T b + c ( 15
)
where the elements of the symmetric, positive-definite matrix A are:

A ef = i∈I -→ W ei • -→ W f i ( 16
)
the rows of vector b are:

b e = i∈I - → v i • -→ W ei ( 17
)
and the scalar c is given by c = i∈I -→ v i 2 . The error E(ρ) on the imposed vessel displacements V is thus minimal when Aρ = b. The optimal set of weights is found in a few iterations using the Conjugate Gradient method, and taking as initial guess ρ = (1, 1, . . . , 1) and the overall deformation U is finally assembled according to Eq. 11.

While the estimated brain model deformation approximates the intraoperatively observed vessel displacements, the reduced number of degrees of freedom in ρ enables a faster computation of the deformation and acts as a smoothing operator strongly limiting the influence of outliers resulting from the elastic registration. As used here, the biomechanical model of the brain can be seen as a "smart interpolator" of the vascular tree deformation.

Planning update and verification

Once the biomechanical model deformation has been computed, the neuronavigator is able to update all the pre-operative data accordingly. The segmented tumor limits defined as surface meshes within the brain FE model are simply deformed by applying the displacements given by the biomechanical model deformation to their vertices.

Pre-operative volumetric images must also be updated in order to reflect the intraoperative changes. Again, the deformed volumetric data is computed by applying to each MR voxel the 3D displacement computed from the deformation of the biomechanical model. The new positions of the resulting voxels are used to rebuild the updated MR volume. This more complex task can be achieved using the 3D texture capabilities of graphic libraries such as OpenGL 9 . The deformed MR volume reconstruction is done in less than one second as 3D texture mapping is a fast procedure carried out within the chipset of most of recent graphics boards.

Among all the updated pre-operative structures, the reference vascular tree plays a crucial role. We describe here how this structure can be used by the surgeon to assess in the region of interest the accuracy of the deformation computed by the system. First the reference vascular tree is deformed in the same way as any anatomical structure defined preoperatively in the MR images, by applying the biomechanical deformation to the segmented centerlines. The deformed structure is called "updated vascular tree". Then, the surgeon performs a freehand Doppler US scan of the region of interest yet this time the US images are not recorded but processed and displayed in real-time. On the processed US images are represented 1) the actual Doppler signal given by the US device and 2) a simulated Doppler signal created using the updated vascular tree. The latter is obtained in the following way. As each US image is localized in space, its position with respect to the brain model can be computed using the transformations described above. The intersection of each US image with the updated vascular tree yields a set of points at the location of the virtually deformed vessels. By comparing the position of these vessels with the actual ones displayed on the same image the surgeon can assess the accuracy of the computed deformation near the lesion and decide whether or not to follow the indications of the system.

Modeling tissue resection

Biomechanical modeling of tissue resection has been widely studied in the literature [START_REF] Vigneron | FEM/XFEM-Based Modeling of Brain Shift, Resection, and Retraction for Image-Guided Surgery[END_REF][START_REF] Jerabovka | Volumetric modeling and interactive cutting of deformable bodies[END_REF][START_REF] Courtecuisse | Gpu-based real-time soft tissue deformation with cutting and haptic feedback[END_REF][START_REF] Cotin | An hybrid elastic model for real-time cutting, deformations, and force feedback for surgery training and simulation[END_REF]. To our knowledge all of the proposed techniques rely on an accurate tracking of the cutting tool (scalpel, aspirator) in order to update in real time the model topology. These techniques are well suited for surgery simulation where the ancillary position is controlled at each time step by the haptic interface. In our experience however it is very difficult to ensure continuous and accurate tracking of the complex surgical gestures involved in tissue resection [START_REF] Bucki | Modélisation Biomécanique des Tissus Mous du Cerveau et Développement d'un Neuronavigateur Permettant la Prise en Compte Per-Opératoire du Brain-Shift[END_REF].

To overcome this issue a tracking-free resection modeling technique is proposed here.

First the surgeon freely performs the resection of the damaged tissues in accordance with the surgical routine. Once tissue removal is complete the shape of the cavity is acquired through intraoperative imaging as discussed in [START_REF] Bucki | Image-based soft-tissues resection tracking for intraoperative brain-shift compensation: general framework and preliminary study[END_REF]. The iterative relaxation procedure described below leads to an equilibrium configuration which complies with both the imposed displacements -stemming from the vascular tree tracking, and the shape of the cavity -scanned intraoperatively. The difficulty of the approach stems from the fact that the shape of the resection cavity is only known in the deformed configuration while the brain mesh has been defined in the pre-operative, underformed configuration.

At each step the 3D shape of the resection superimposed on the model virtually separates the undeformed brain mesh in two regions: a region where the elements are outside the resection and a region where elements are inside the resection. Prior to relaxation the nodal forces are assembled so that the elements inside the resected area do not affect the tissues modeled by the elements found outside the resection. Interactions between elements lying fully inside or outside the resection are computed normally, the action-reaction principle is only broken at the resection boundary.

The relaxation of the two sub-meshes defined above are computed alternatively. After each iteration some elements may transit through the resection boundary becoming virtually resected or, to the contrary joining the non-resected region. The iterations stop once the model has reached an equilibrium in the non resected region and no element transition have been observed.

Finally the undeformed configuration of the resection is computed by attaching the cavity surface to the mesh elements and then replacing each mesh node in its initial, preoperative position. The undeformed configuration of the resection cavity can be overlaid on the planned resection contours which makes it possible for the surgeon to assess the extent and location of over-or under-resected areas. Fig. 4 illustrates this technique on a phantom model of the brain and a resection surface extracted from localized 2D US image analysis.

The technique has been tested on a dataset -vessel displacements and resection surface -acquired on a silicone phantom with imaging conditions close to the ones found in an operating room. The relaxation algorithm converges in a couple of minutes yielding an undeformed configuration of the resection cavity with a submillimetric mean accuracy (maximal errors are the order of the millimeter) [START_REF] Bucki | Image-based soft-tissues resection tracking for intraoperative brain-shift compensation: general framework and preliminary study[END_REF].

A clinical case

This section illustrates our neuronavigation technique on a clinical case. Section 5.1 describes the pre-operative data acquisition and analysis. The brain-shift measurements performed during the intervention are discussed in §5.2. Finally, section 5.3 gives the details of the postoperative measurements revealing the errors committed by the system on the brain-shift estimation. The neuronavigator prototype presented here has not been certified for clinical use and, although all the computations were performed during the intervention, the output was hidden to the surgeon so as to not affect his decisions. The evaluation of the system discussed below was performed on a 46 year old female patient with a large meningioma (approx. 125 cm 3 ) in the right temporal lobe. The patient presented clinical signs of intracranial hypertension probably due to the tumor growth. Before the intervention an informed consent was obtained and signed by the patient.

Pre-operative phase

Two MR volumes of the patient's head were acquired before the intervention. A 512 × 512 × 120 T1 MR volume, with 0.39 × 0.39 × 1.5mm voxel resolution, used for the routine neuronavigation was recorded after 7 skin fiducials were placed on the patient's face. We will refer to this image as the "Surgiscope" volume.

A second 512 × 512 × 120 volume with a 0.39 × 0.39 × 0.7mm voxel resolution was acquired after the patient was injected with Gadolinium for better visualization of distal brain vessels. We will refer to this image as 3DToF MRA or simply MRA. The covered volume was approximately half of the Surgiscope volume and was focused on the region around the tumor. Fig. 5-a presents an axial slice from the Surgiscope MR volume showing the hyperintense lesion. On the corresponding MRA slice, shown in Fig. 5-b, the brain vascular tree and more specifically the tumor vascularisation appear as bright spots.

Patient skin surface was extracted from both the Surgiscope and MRA volumes using a Marching Cubes procedure [START_REF] Lorensen | Marching cubes: a high resolution 3-d surface construction algorithm[END_REF] and the corresponding surfaces were rigidly registered by an Iterative Closest Point algorithm [START_REF] Besl | A method for registration of 3d shapes[END_REF] in order to align the volumes in a common image referential. The tumor was segmented in the Surgiscope volume with the 3D Slicer software using the simple region growing technique after a seed was manually placed in the hyperintense lesion region. A smooth surface model of the tumor was reconstructed from the resulting binary volume by elastically deforming and fitting the template shape of a bounding ellipsoid on the points cloud formed by the boundary voxels, i.e. the segmented tumor voxels in a 26-neighborhood of a non-tumor voxel. The lesion volume was estimated to be 125 cm 3 .

The 3DToF MRA volume segmentation, on the other hand, produced a set of 2612 vessel centerlines defining the patient's vascular tree pre-operative configuration. Fig. 5-c shows a close-up of the reconstructed tumor volume, shown as a transparent surface, aligned with the points cloud of the reference vascular tree. Finally a patient specific FE mesh was generated following the procedure described in §3.2, after the surgical planning was established and entered on the Surgiscope console by the surgeon.

The generic FE mesh was registered onto the patient's cortex within 20 seconds producing the patient specific model shown in Fig. 7-b. The position of the craniotomy defined by the surgeon was used to identify the cortical nodes lying inside the surgical opening on the free surface of the brain. As discussed in §3.2, frictionless brain-skull sliding was restricted to the tangent plane at all cortex nodes on the lower half of the brain. Finally the brain stem nodes were fixed.

The set of segmented brain vessels was localized within 713 elements of the FE mesh. The deformation precomputations were carried out only on the 482 elements located in the hemisphere undergoing surgery as the limited US penetration depth was unlikely to reveal the vascularisation on the opposite side. The 4 × 3 × 482 = 5784 precomputed deformation vectors, one for each spatial dimension and each node of each involved tetrahedron, were estimated in about 8 hours (1 minute processing per-element), and the 170 megabytes of resulting data (30 kilobytes per vector, stored as 4 bytes double precision floating point values) were saved on a CD-ROM to be retrieved at system start-up on the day of the intervention.

Intraoperative phase

Once the patient was installed in the OR, the positions of the 7 skin fiducials were acquired using a calibrated pointer and converted into the patient's head referential. The physical marker positions were registered with the ones identified in the Surgiscope MR volume. The accuracy of the Arun registration of the fiducials was: mean error=0.8 mm, max error=2.1 mm and RMS=0.4 mm.

After craniotomy and dura opening a brain sagging could be observed and a series of US images of the brain were acquired in order to capture the brain-shift. To this end, the craniotomy cavity was filled with sterile physiological fluid and a continuous 4 minute-long freehand Doppler US sweep of the region of interest was performed by the surgeon. The 715 localized frames were segmented within 15 seconds yielding a set of 1920 vessel centerline points.

The elastic registration and filtering required approximately 5 seconds and produced a set of 156 vessel displacements. The statistics on the displacement field were: mean displacement=6.0 mm, maximal displacement=9.8 mm, RMS=6.2 mm. Figure 6 illustrates the brain vessels registration procedure described in §4.4. The biomechanical model was deformed according to the estimated brain vessels displacements field as described in §4.5. The number of elements involved in the deformation computation was N E = 65. The optimal weight values ρ and the corresponding FE model deformation were computed in less than 2 seconds. The error on the 156 imposed vessel displacements was: mean error=1.5 mm, max. error=7.3 mm and RMS=2.0 mm. The resulting volumetric deformation was applied to the pre-operatively acquired Surgiscope MR volume, segmented tumor contours and reference vascular tree in 150 milliseconds. Fig. 7-a shows a sample axial MR slice taken from the intraoperatively computed MR volume after the brain-shift correction has been applied. Fig. 7-b shows the intraoperatively deformed FE mesh and tumor surface model. Fig. 7c shows a close-up on the region where the tumor margin exhibits the maximal 6mm shift. The boundaries of the preoperative and deformed tumor margins are superimposed using transparency which clearly shows the margin displacement produced by the brain shift. The dark dots represent the manually segmented US points that lie on the deformed tumor margin and which were used to assess the accuracy of the brain-shift compensation. The image shows a much better adequacy for the corrected tumor margin (rightmost) than for the preoperatively segmented margin (leftmost).

Postoperative deformation control

Immediately after the deformation computation a set of 503 US B-mode images has been acquired to measure the error in the estimation of brain deformation.

After the intervention an operator manually segmented the images where the deeper margin of the lesion was clearly visible by placing points at the interface between the meningioma and the cortex. A set of 1246 points localized on the deformed tumor margin was obtained by converting the segmented points into the patient's head reference frame using US probe calibration information along with the position of each recorded US image.

The point-to-surface distance between the resulting points cloud and the tumor limits was measured. Table 1 gives the error measurements on tumor localization before and after brain-shift compensation. The intraoperative compensation for the brain-shift deformation resulted in a significant gain of accuracy as 75% of the deformation was captured by the system.

Mean

Max RMS Before compensation 3.5 mm 7.6 mm 3.7 mm After compensation 0.9 mm 1.7 mm 1.0 mm Table 1. Error on tumor margin localization before and after brain-shift compensation.

Conclusion and Perspectives

A prototype of a neuronavigation system featuring fast and fully automatic intraoperative brain-shift compensation has been presented. It relies on the intraoperative tracking of the brain vascular tree deformation through localized 2D Doppler US imaging and patient-specific FE modeling of the brain soft tissues. A partial brain deformation field is inferred by elastic registration of the undeformed pre-operative and deformed intraoperative vascular tree configurations. The deformation is further extended to the whole organ using a patient specific FE biomechanical model of the brain. The system is very reactive as it is able to deliver a deformation estimation along with updated pre-operative data in about 10 seconds which ensures that the surgeon is provided with relevant information about the configuration of the patient's anatomy. However this reactivity comes at the cost of a heavy precomputation that must be run before the intervention, as soon as medical images of the patient's head have been acquired.

During the evaluation of the system in clinical conditions on a patient presenting a large meningioma the system was able to compensate for 75% of the brain-shift. Although a wider study is necessary to assess the system performance, this preliminary work suggests adequacy of the developed techniques to address this issue. From a practical point of view, intraoperative US imaging was performed without raising any major technical issues or additional risks for the patient, and the use of sterile physiological fluid as coupling medium proved to be effective. Unlike interventional MRI Doppler US is a low-cost imaging modality readily available in most neurosurgical operating rooms, and furthermore, Doppler US acquisitions can easily be repeated a number of times during the surgical procedure in order to capture the dynamics of brain-shift. Finally, thanks to the color encoding of the flow information in the US images, the segmentation of the acquired images could be performed automatically within a couple of seconds.

As for the patient specific biomechanical modeling, an assumption was made that the quantity and quality of intraoperative information on tissue displacements retrieved from Doppler US imaging could compensate for the limited knowledge of the tissues rheology. Indeed, as the newly formed brain vessels are more specifically localized near the tumor margin, the displacement boundary conditions resulting from the vascular tree tracking strongly constrain the deformation computation. As a consequence, the FE model plays a less central role here than in systems relying solely on cortical surface tracking. This hypothesis as well as the need for fast and robust deformation computation led to the choice of the linear (small strain and small displacements) framework to model the brain soft tissues.

It is important to note that the clinical case presented here made all measurements possible as the treated meningioma margin was clearly visible in the control B-mode US images. However, the expected application field of the system are situations where tumor limits are hard to identify by palpation or US imaging. Unfortunately, in these cases, system validation is hard to perform without relying on intraoperative MR imaging precisely because this modality is the only one that allows accurate tumor margin tracking.

Concerning patient pre-operative imaging, 3DToF MRA is a delicate technique requiring experimentation and fine tuning of the acquisition parameters to reach optimal imaging quality. Its main limitations are due to signal loss near flow irregularities, vessels obstruction by superposition with fat structures, spin saturation induced by slow flow in distal arterial branches and patient motion during the usually long acquisition time [START_REF] Wilms | Magnetic resonance angiography of the intracranial vessels[END_REF]. Furthermore, high dose contrast agent administration (0.01 -0.02 mmol/kg) may result in excessive enhancement of the soft tissues and the venous system in the images especially around the Circle of Willis, and may thus obscure the arteries [START_REF] Jung | Contrast-enhanced MR angiography for the diagnosis of intracranial vascular disease: optimal dose of gadopentetate dimeglumine[END_REF]. Although alternate methods, such as the selective water excitation technique, have been proposed in the literature to enhance vessels visualization by eliminating the stationary background from the acquired images [START_REF] Őzsarlak | 3D time-of-flight MR angiography of the intracranial vessels: optimization of the technique with water excitation, parallel acquisition, eight-channel phased-array head coil and low-dose contrast administration[END_REF], none of them was put in practice during this preliminary study and only routine 3DToF MRA protocols were used.

The identification of the reference vascular tree is a crucial step in our procedure as the brain-shift estimation relies on the tracking of the vascular tree deformation. The semi-automatic procedure used here is not optimal and could be enhanced using region growing techniques as suggested in [START_REF] Descoteaux | A multi-scale geometric flow for segmenting vasculature in mri: Theory and validation[END_REF]. As for the intraoperative Doppler US acquisition the 3Hz framerate is due to the performance of the video analogic-to-digital conversion board. This results in a lengthy recording time which could be greatly shortened by using 3D ultrasound probes and software featuring signal post-processing that allows reconstruction of continuous volumetric US images, as opposed to the pseudo-volumes used here. Furthermore a single-shot 3D ultrasound acquisition performed with a steady probe, as opposed to our freehand sweep, would strongly reduce the number of Doppler artifacts created by unwanted tissue displacements [START_REF] Ji | Coregistered volumetric true 3d ultrasonography in image-guided neurosurgery[END_REF]. Most importantly, wide US probe movements are likely to affect the brain tissue configuration and impair the accuracy of the computation. For the above mentioned reasons switching to 3D ultrasound is the next step that has to be taken in order to improve the neuronavigation system described here. In this regard 3D ultrasound devices can reasonably be expected to be found in the operating theaters in the near future.

Compared to the technique developed by Reinertsen et al., the approach proposed here relies on biomechanical modeling and is capable of taking into account brain-skull contact boundary conditions, heterogeneous tissue behavior and complex interactions such as tissue retraction or resection which seem much more complex to implement using a TPS formulation. Tumor resection is the cause of most intraoperative changes occurring during surgery and therefore, its modeling is the upcoming challenge of model-based neuronavigation. The original approach for intraoperative tissue resection tracking and FE modeling described above will soon be evaluated and integrated in the neuronavigator along with the broadening of the modeling framework to account for both mechanical and geometrical non-linearities.

Finally, from a more theoretical point of view, the decision as to when to perform a brain-shift update rests with the surgeon and thus is limited by the physician experience and subjective perception of the brain deformation. The planning update could alternately be triggered by objective criteria such as the cortex movements measured by a laser range scanner or a pointer but this approach ignores subcortical changes and seems delicate to put in practice once tissue retraction or resection has been performed. A more suitable solution would consist in leaving a small transducer in the cranial opening to continuously scan the region of interest and monitor the vessel displacements. Coupled with the neuronavigator capabilities, this "brain-shift monitoring" system would be able to update the surgical planning every 20 to 30 seconds. Although a bit early, constant progress in 2D and 3D endosonography [START_REF] Lorenzo | Tridimensional (3d) endoscopic ultrasonography[END_REF] [66] [START_REF] Saftoiu | Tridimensional (3d) endoscopic ultrasound a pictorial review[END_REF] brings this proposal closer to feasibility.

Fig. 1 .

 1 Fig. 1. Neuronavigator workflow during the pre-and intraoperative phases.

Fig. 2 .

 2 Fig. 2. Neuronavigator elements layout in the operating room.

Fig. 3 .

 3 Fig. 3. Elastic registration filtering. The source points are represented by gray disks and destination points as black disks. In (a), the computed elastic deformation R brings all source points near D within d max (disks radius in the figure) except for the white disk, which displacement (dashed arrow) does not comply with the forward filtering criterion. The corresponding source point is thus removed from the set. In (b), the inverse registration function R -1 is applied to all destination points. Only the full black disks reach S f within accuracy d max . The resulting displacement vectors set V, stemming from D b and pointing to a subset of S f , is represented in (c).

Fig. 4 .

 4 Fig. 4. Computing the preoperative shape of an intraoperatively identified resection cavity. (a) Cutout view of the preoperative FE brain model overlaid with the intraoperative resection cavity (black region). (b) Brain model shifted downwards due to gravity; the final steady state accounts for the absence of tissue in the resected region. (c) Brain model and resection cavity in the preoperative frame (black region) overlaid with the cavity in the intraoperative frame (grey region).

Fig. 5 .

 5 Fig. 5. (a) Axial slice taken from the Surgiscope volume. (b) Equivalent slice in the 3DToF MRA image. (c) Reconstructed 3D surface of the lesion and reference vascular tree segmented from the 3DToF MRA.

Fig. 6 .

 6 Fig. 6. Intraoperative brain vessels registration and displacements field computation. (a) Doppler vessels (gray) and pre-operative vascular tree (black) before elastic registration. (b) After registration and filtering. (c) Resulting displacements field.

Fig. 7 .

 7 Fig. 7. Updated pre-operative data and tumor margin shift. The arrows indicate the upward direction. (a) A sample axial slice from the intraoperatively recomputed MR volume. (b) Patient specific FE brain model, segmented tumor surface and estimated vessel displacements field. (c) Tumor margin before (left) and after (right) brain-shift deformation.

ISIS, Saint Martin d'Héres, France.

www.slicer.org

Northern Digital Inc., Canada.

Praxim, La Tronche, France.

The Imaging Source Europe GmbH, Germany.

http://www.opengl.org

Acknowledgements

The authors wish to thank Dr. Pierre Lavagne and Dr. Gilles Francony from the surgical reanimation unit (URC) at the Hospital Michallon for their advice and help during this study.