M Lau

X Faure

F Zara

F Jaillet

J-M Moreau

Posters and Demos

Keywords: (according to ACM CCS): I, 3, 5 [Computer Graphics]: Computational Geometry, Object Modeling-Physically based modeling I, 6, 8 [Simulation And Modeling]: Types of Simulation-Parallel

Figure 1: (1) Rendered beam for Hooke's and Saint Venant-Kirchhoff's material, and initial state (from left to right) -(2-4) Deformation simulation of a rabbit (initial 3D mesh courtesy of L. Stanculescu).

Introduction

The Tensor-Mass (TM) approach is a good candidate to handle deformable objects within the context of physicallybased simulations in complex interactive scenes. The deformation forces are derived from the FEM mechanical formulation, but are next computed locally and iteratively for each discretized element. Besides, this allows to handle more directly and easily topological changes (cutting) and external interactions (collisions, tools), in the same way as Mass-Spring Systems.

Hence, several formulations have been proposed to account for various mechanical behaviors: the linear Hookean model [START_REF] Cotin | A hybrid elastic model for real-time cutting, deformations, and force feedback for surgery training and simulation[END_REF]; the non-linear geometrical model based on Saint Venant-Kirchhoff's elasticity model [PDA00], or anisotropic material [Pic03]; or an extension with pre-computation for non-linear visco-elastic deformations [SDR ⇤ 05].

But, as far as we know, the problem of the parallel implementation of the TM model has not yet been addressed. In this paper, a GPU implementation is presented, considering both linear and non-linear mechanical behaviors and using an implicit integration scheme to ensure the unconditional stability of our simulation.

Simulation of the deformation of an object

The TM approach is based on the domain's discretization into several elements, and the mechanical equations are next solved locally, involving the following main steps for each element:

• Discretization of the displacement U E with the definition of interpolation functions L according to the chosen type of elements (hexahedron, tetrahedron, etc.); • Computation of the strain-tensor according to the chosen elasticity model amongst the most employed in interactive simulations, namely a) Hooke's or b) Saint Venant-Kirchhoff's elasticity models -within the element, at X:

a) e l (X) = 1 2 (rU T (X) + rU(X)) b) e nl (X) = 1 2 (rU T (X) + rU(X) + rU T (X) rU(X))
• Computation of the deformation energy W E , according to the displacement of the element's node;

c The Eurographics Association 2012.

X. Faure & F. Zara & F. Jaillet & J-M. Moreau / Implicit Tensor-Mass solver on the GPU

• Computation of the elasticity force F E by the evaluation of the derivative of the deformation energy W E .

Force computation on the GPU

The parallel computation of the forces on the GPU is divided into 2 tasks, involving a set of kernels (Alg. 1):

1. First, we compute and store the forces applied on each node of a given element. This computation (kernel1) is made for each element of the domain. 2. Next, we sum these partial forces to obtain the total forces applied on each node, involved by the different elements of the domain. This computation (kernel2) is made for each node of the domain.

Algorithm 1 Force computation on the GPU 1: {N : number of elements; m : total number of nodes} 2: {n : number of nodes per element} 3: {Nn : max number of neighbor elements for a node} 4: // Task 1-Computation of partial forces 5: for e = 0 to N 1 do 6:

// Execution of N kernel1 7:

for v = 0 to n 1 do 4. Implicit integration method on the GPU Then, the dynamical equations are derived from Newton's laws to compute the object acceleration. An implicit integration scheme, requiring computation of the derivatives of forces, is then used to obtain deformation and displacement. The parallel algorithm for this is similar to Alg. 1with the only difference in the use of DForce() instead of Force(), which depends on the same parameters but also on the time step h and the current nodes velocity, to directly compute h 2 ∂F ∂U V (t) or h ∂F ∂U DV , and avoid the storage of the sparse Jacobian matrix ∂F ∂U . Next, the linear system is solved on the GPU with the Conjugate Gradient method, following [START_REF] Allard | Implicit FEM Solver on GPU for Interactive Deformation Simulation[END_REF].

Results and performances

We compare running times on the CPU and GPU, using Linux Ubuntu11.04. The CPU is an Intel R Xeon R 4 cores @3.07 GHz. The GPU is a GeForce GTX 560, 2047 MB, 56 cores @1.620 GHz. Fig. 2 presents the speedup (ratio between GPU and CPU execution time) obtained for triangles and tetrahedra for both the linear and non-linear TM model (Fig. 1). Moreover, similar results are presented for SOFA's FEM implementation, i.e. the corotational FEM model for tetrahedra [NP05], in its GPU version according to [START_REF] Allard | Implicit FEM Solver on GPU for Interactive Deformation Simulation[END_REF]. A speedup of 25.5 is reached for SOFA's FEM implementation and of 29.5 for our TM, for a beam composed of 307,200 elements.

Conclusion and perspectives

In this paper, we presented an original implementation of the TM model on the GPU that considerably speeds up the simulation times. Comparisons between running times on the CPU and the GPU suggest that the parallel implementation of the model becomes interesting for increasingly complex computations.

  [e][v]][index[e][v]]=Force(...); 9: end for 10: end for 11: // Task 2 -Sum of partial forces 12: for i = 0 to m do 13: // Execution of m kernel2 14: for j = 0 to Nn 1 do 15: TotalForce[i] += PartialForce[i][j]; 16: end for 17: end for Moreover, specific data structures are defined for a discretization of the domain into N elements involving m nodes: • index (of size N ⇥ n) stores the relationship between local and global indexation for each node of each element. • PartialForce (of size Nn ⇥ m ⇥ 3) enables the storage of 3D coordinates of partial forces for each node considering its global indexation. • TotalForce (of size m ⇥ 3) stores the sum of the partial forces for each node considering its global indexation. • ForceIndex (of size N ⇥ n) stores the index of data structure PartialForce.

Figure 2 :

 2 Figure 2: Performances of the TM parallelization.

c The Eurographics Association 2012.


B Denis De Senneville

B Quesson

C T W Moonen
email: moonen@imf.u-bordeaux2.fr

Magnetic Resonance Temperature Imaging

Continuous, real-time 3D temperature mapping during a hyperthermic procedure may provide enhanced safety by visualizing temperature maps in and around the treated region, improved efficiency by adapting local energy deposition with feedback coupling algorithms, and therapy endpoints based on the accumulated thermal dose. Noninvasive mapping of temperature changes can be achieved with MRI, and may be based on temperature dependent MRI parameters. The excellent linearity of the temperature dependency of the proton resonance frequency (PRF) and its near-independence with respect to tissue type make the PRF-based methods the preferred choice for many applications, in particular at mid to high field strength (≥ 0.5 T). The PRF methods employ RF-spoiled gradient echo imaging methods, and incorporate fat suppression techniques for most organs. A standard deviation of less than 1 o C, for a temporal resolution below 1 s and a spatial resolution of about 2 mm, is feasible for immobile tissues. Special attention is paid to methods for reduction of artifacts in MR temperature mapping caused by intra-scan and inter-scan motion, and motion and temperature-induced susceptibility effects in mobile tissues. Real-time image processing and visualization techniques, together with accelerated MRI acquisition techniques, are described because of their primary importance for real-time, image guided, therapy guidance.

INTRODUCTION

Heat conduction through diffusion and perfusion processes vary locally as a function of tissue architecture, tissue composition, physiological parameters and temperature. Energy deposition at or near the target during thermal therapies depends on the efficiency of absorption of the transmitted energy employing electromagnetic (optical, radiofrequency, microwave) or acoustic waves, and is therefore also a function of tissue composition. Tissue coagulation during the procedure may significantly modify local heat conduction as well as energy absorption. Continuous thermometry appears therefore helpful for assuring adequate thermal treatment in the target region, and avoiding damage to healthy regions. For RF, laser or microwave heating, a needle is inserted in or near the target tissue, and a thermometer is often incorporated in the needle. This allows thermometry at a single point, but does not provide further information on the spatial distribution of temperature for the treatment of a larger volume around the heat source. In case of FUS, the transducer is far from the focal point, generally even outside the body. An inserted thermometer interferes with the ultrasound waves with the wavelengths that are typically used for FUS hyperthermia.

Therefore, classical thermometry sensors offer limited value for thermotherapy.

Of the different imaging modalities, MRI appears the ideal tool for temperature mapping. A particular advantage of MRI for guiding thermal procedures is that MRI not only allows temperature mapping but it can be used as well for target definition, and may provide an early evaluation of the therapeutic efficacy. Whereas temperature effects on MR measured parameters (resonance frequency, relaxation, diffusion) had been well described previously, the first report of temperature mapping by MRI appeared in 1983 [START_REF] Parker | Temperature distribution measurements in two-dimensional NMR imaging[END_REF]. The method was based on the longitudinal relaxation time (T 1 ). Since then, alternative MR temperature imaging methods have been proposed, based on diffusion coefficient (D) [START_REF] Bihan | Temperature mapping with MR imaging of molecular diffusion: application to hyperthermia[END_REF], proton resonance frequency (PRF) (3) of tissue water and other temperature-dependent parameters. In this paper, the different methods for MR temperature mapping are reviewed. Special attention is paid to sensitivity, and linearity of the effects with temperature, their dependence on coagulation, dependence on field strength, and the question whether absolute or relative temperature is measured. In addition, pulse sequences, speed, potential artifacts, and motion sensitivity are discussed for each method, in particular with respect to thermotherapy procedures. Recently, progress in temperature mapping has been reported using ultrasound imaging (4) but these alternative methods still require significant refinement before they can be used for guidance of hyperthermic therapies.

1.

Physical basis of temperature MRI methods

Temperature MRI based on the water proton resonance frequency

The temperature dependence of the water proton resonance frequency (PRF) was first investigated by Hindman for the study of intermolecular forces and hydrogen-bond formation between water molecules [START_REF] Hindman | Proton resonance shift of water in the gas and liquid states[END_REF]. MR temperature mapping on the basis of the PRF was first proposed by Ishihara et al. [START_REF] Ishihara | A precise and fast temperature mapping method using water proton chemical shift[END_REF] and further developed by de Poorter [START_REF] Poorter | Noninvasive MRI thermometry with the proton resonance frequency method: study of susceptibility effects[END_REF][START_REF] De Poorter | Noninvasive MRI thermometry with the proton resonance frequency (PRF) method: in vivo results in human muscle[END_REF]. It is currently the most widely used MR temperature mapping method at mid and high field (0.5T and above).

The theory is briefly explained below.

The local magnetic field B nuc as observed by the spins is a function of the main magnetic field B 0 and the temperature T dependent chemical shift field σ:

0 )) ( 1 ( ) ( B T T B nuc σ + = [1]
The chemical shift field (in ppm) is the sum of temperature independent contributions, for example those originating from B 0 field inhomogeneities, represented by σ 0 , and a temperature dependent contribution, σ T (T):

) ( ) ( 0 T T T σ σ σ + = . [2]
The chemical shift field can be calculated from the phase information in RF-spoiled gradient

echo images: 0 ) ( ) ( B T T T E γσ = Φ [3]
where φ is the image phase, γ is the gyromagnetic ratio of the observed nucleus (42. 58 10 6 Hz T -1 for protons), and TE is the echo time. In order to measure temperature dependent changes in chemical shift, the term σ 0 must be eliminated, which is typically accomplished by subtraction of the field distribution measured at a given reference temperature from the field distribution measured at temperature T, leading to:

0 ) ( ) ( B T T T T T T E ref ref αγ Φ - Φ = - = ∆ [4]
where α is the temperature dependency of the water chemical shift in ppm⋅°C -1 .

In principle, any gradient echo method can be used for PRF-based MR thermometry, so long as contributions from stimulated echoes can be neglected [START_REF] Chung | Signal formation in echo-shifted sequences[END_REF][START_REF] De Zwart | Fast lipid-suppressed MR temperature mapping with echo-shifted gradient-echo imaging and spectral-spatial excitation[END_REF]. RF spoiling of fast gradient echoes is thus necessary when flip angles close to the Ernst angle are used for optimal SNR for short TR. Spin-echo sequences cannot be used since the temperature induced phase contribution will be refocused.

The phase difference between reference data and data acquired during heating increases linearly with the echo time of the experiment, whereas the image signal-to-noise ratio decreases exponentially as a function of TE. Therefore the optimal echo time of the thermometry sequence is equal to the T 2 *-value [START_REF] Chung | Optimization of spoiled gradient-echo phase imaging for in vivo localization of a focused ultrasound beam[END_REF][START_REF] De Zwart | Fast magnetic-resonance temperature imaging[END_REF]. Resulting SNR can be expressed as:

* 2 T E T E e cT SNR - = [5]
where c is a constant, depending on spin density, flip angle, TR, as well as instrumental parameters (field strength, RF coil).

The temperature dependence of the water proton resonance frequency is approximately 0.01 ppm⋅°C -1 . A very important advantage of the PRF method is its nearindependence of tissue composition [START_REF] Peters | Ex vivo tissue-type independence in protonresonance frequency shift MR thermometry[END_REF]. The small remaining tissue dependency may be related to ion concentration [START_REF] Peters | Magnetic resonance thermometry for predicting thermal damage: an application of interstitial laser coagulation in an in vivo canine prostate model[END_REF]. The proportionality constant of PRF with temperature appears at first sight rather low. Comparisons between T 1 , D, and PRF based methods have resulted in higher precision for PRF methods [START_REF] Wlodarczyk | Threedimensional monitoring of small temperature changes for therapeutic hyperthermia using MR[END_REF]. An example of PRF-based temperature imaging in vivo is shown in Figure 1 (from ( 16)).

It has been assumed in PRF methods that the macroscopic magnetic susceptibility is independent of temperature, as shown in Equation [START_REF] Bihan | Temperature mapping with MR imaging of molecular diffusion: application to hyperthermia[END_REF]. However, it has been demonstrated that temperature dependent changes in magnetic susceptibility are not negligible [START_REF] Poorter | Noninvasive MRI thermometry with the proton resonance frequency method: study of susceptibility effects[END_REF][START_REF] Peters | Heat-source orientation and geometry dependence in proton-resonance frequency shift magnetic resonance thermometry[END_REF][START_REF] Young | An evaluation of the effects of susceptibility changes on the water chemical shift method of temperature measurement in human peripheral muscle[END_REF]. It has been suggested that such effects can be corrected in real-time [START_REF] Salomir | A fast calculation method for magnetic field inhomogeneity due to an arbitrary distribution of bulk susceptibility[END_REF]. The presence of lipids is a potential source of artifacts since the PRFs of lipid hydrogens is independent of temperature. Lipids can be conveniently suppressed in gradient echo imaging by frequencyselective slice excitation [START_REF] De Zwart | Fast lipid-suppressed MR temperature mapping with echo-shifted gradient-echo imaging and spectral-spatial excitation[END_REF] or alternative methods. PRF-based methods are also sensitive to microscopic susceptibility effects, e.g. related to changes in the concentration of deoxyhemoglobin during the therapy, similar to the BOLD effect [START_REF] Ogawa | Brain magnetic resonance imaging with contrast dependent on blood oxygenation[END_REF]. Also, strong temperature gradients within a single voxel will lead to a large phase dispersion and hence signal decrease.

Apart from SNR advantages at high field, the sensitivity of the PRF-method increases linearly with the magnetic field strength. However, the exact relationship between sensitivity of PRF-based thermometry and magnetic field strength is not straightforward since gradient echo images are acquired with long TE. The T 2 * -value of the tissue generally decreases with increasing field strength because of macroscopic and microscopic susceptibility effects. Since T 2 * contains T 2 and susceptibility contributions, the relation of T 2 * with field strength depends on the relative contributions of inhomogeneity and T 2 to the observed T 2 * . If inhomogeneities are dominating, the T 2 * will decrease approximately linearly with increasing field strength, whereas a T 2 * -value close to the T 2 -value of the tissue will hardly depend on field strength (assuming T 2 is rather independent of field strength). Note that macroscopic susceptibility effects can be decreased by improving the spatial resolution. An additional small magnetic field effect originates from the increase in T 1 relaxation time with increasing field resulting in a lower Ernst-angle and in a decreased image SNR.

Non-invasive PRF-based temperature imaging only accounts for relative temperature changes, by comparing a reference image to images acquired under identical experimental conditions. As a consequence, absolute temperature cannot be evaluated. An alternative approach to calculate absolute temperature is MR spectroscopic imaging to measure water resonance shift. The principles of the method are similar to those of PRF with the difference being the acquisition method. Small frequency drifts of the instrument would induce substantial errors on the PRF calculated temperature values whereas proton spectroscopic imaging could include a reference frequency (e.g. lipids [START_REF] Kuroda | Temperature mapping using the water proton chemical shift: a chemical shift selective phase mapping method[END_REF], N-acetyl-aspartate in brain ( 22))

and thus be relatively immune to field drifts.

MR temperature mapping based on other physical mechanisms 1.2.1 Temperature MRI based on the T 1 relaxation time of water protons

The variation in T 1 with temperature can be described by the following relationship (see complete references in [START_REF] Quesson | Magnetic resonance temperature imaging for guidance of thermotherapy[END_REF]).

kT T a E e T T ) 1 ( 1 1
) ( -∞ = [START_REF] Ishihara | A precise and fast temperature mapping method using water proton chemical shift[END_REF] where E a (T 1 ) is the activation energy of the relaxation process, k is the Boltzmann constant and T is absolute temperature. Exchange processes between mobile bulk water and relatively immobilized water at surfaces of proteins and membranes play an important role. Within a small temperature range T 1 is linearly dependent on temperature. However, both the T 1 and its temperature dependence differ between tissues. Non-linear effects have been observed in particular when coagulation occurs [START_REF] Fried | Dynamic T1-weighted magnetic resonance imaging of interstitial laser photocoagulation in the liver: observations on in vivo temperature sensitivity[END_REF][START_REF] Graham | Time and temperature dependence of MR parameters during thermal coagulation of ex vivo rabbit muscle[END_REF][START_REF] Graham | Analysis of changes in MR properties of tissues after heat treatment[END_REF]. This is not surprising since exchange processes and effective tumbling rates of water may be influenced drastically. Typical changes of T 1 with temperature are of the order of 1%/ o C (23). T 1 generally increases from low (<0.5T) to high field (>1T) [START_REF] Bottomley | A review of normal tissue hydrogen NMR relaxation times and relaxation mechanisms from 1-100 MHz: dependence on tissue type, NMR frequency, temperature, species, excision, and age[END_REF], and T 1 contrast between tissues diminishes. Temperature dependency of T 1 also becomes smaller at higher field strength. Apart from SNR advantages at high field (SNR increases approximately linearly with field strength), temperature mapping based on T 1 is more sensitive at relatively low field strength. T 1 methods generally require more time than PRF based methods. A qualitative overview of the specifications of temperature MRI based on T 1 is given in REF [START_REF] Quesson | Magnetic resonance temperature imaging for guidance of thermotherapy[END_REF].

Temperature MRI based on the molecular diffusion constant of water

The thermal Brownian motion of an ensemble of molecules is characterized by the diffusion constant D. The relationship between temperature and the diffusion constant is

exponential kT D a E e D ) ( - ≈ [7]
where E a (D) is the activation energy of the molecular diffusion of water, k is the Boltzmann constant and T is absolute temperature. The dependence of D on temperature can be described .

[8]

As compared to T 1 , the activation energy for diffusion of water is substantially higher and the temperature dependency of changes in D amounts to about 2%/°C. Diffusion constants can be measured by MR on the basis of the signal attenuation observed in the presence of a pair of strong mutually counteracting gradients [START_REF] Delannoy | Noninvasive temperature imaging using diffusion MRI[END_REF].

Despite the good sensitivity of the D method (23), practical problems limit the routine use of the method in vivo. The mobility of water in tissues depends on barriers such as cellular structures, proteins and membranes. This has two major consequences for temperature MRI based on D. Firstly, the dependence of diffusion on temperature is non-linear since the permeability of such barriers is temperature dependent. Coagulation processes lead to large changes in diffusion constant since diffusion pathways are modified. Non-lethal physiological effects can also lead to large changes in D, for example those observed immediately following ischemia in the brain (D is lowered about 40%) when such changes are still reversible [START_REF] Moseley | Early detection of cerebral ischemia in cats: comparison of diffusion-and T2-weighted MRI and spectroscopy[END_REF]. The latter effect has been attributed in part to cell swelling. Similar shifts of water between extracellular spaces and intracellular compartments may also influence D during hyperthermia procedures. Secondly, the mobility of water in tissue is direction-dependent because of the anisotropic nature of the barriers, e.g. in muscle fibers, necessitating full characterization of the diffusion tensor.

It should be noted that diffusion-based temperature MRI methods are based on diffusional displacements of about 10 micrometers via the use of very strong gradient pulses.

The methods are therefore also sensitive to macroscopic motion. Apart from SNR advantages at high field and small effects of relaxation times on speed of the imaging methods, temperature mapping based on D is independent of field strength.

Temperature MRI based on changes in bulk magnetization and magnetization transfer rates

The use of some alternative thermometry methods have been suggested, e.g. spin density or magnetization transfer (see [START_REF] Quesson | Magnetic resonance temperature imaging for guidance of thermotherapy[END_REF], and references therein). MR signal strength depends linearly on the magnetization M 0 . Since M 0 depends on the Boltzmann thermal equilibrium, it is possible to evaluate temperature changes based on M 0 measurements.

However, the temperature dependency is rather small. In addition, when using fast MRI, modifications in M 0 are difficult to separate from effects due to changes in relaxation times.

Magnetization transfer depends on chemical exchange of bulk water and labile chemical groups at the surface of proteins and membranes. Such exchange process is temperaturedependent and may be used for MR thermometry. However, the sensitivity of these methods appears limited and strongly tissue-dependent, and no full analysis is therefore presented in this paper.

Temperature MRI based on temperature-sensitive contrast agents

The temperature sensitivity of some compounds is much higher than that of water. A recent report noted a temperature dependence of more than 1 ppm/ o C [START_REF] Aime | A new ytterbium chelate as contrast agent in chemical shift imaging and temperature sensitive probe for MR spectroscopy[END_REF]. If such agents can be introduced without toxic side effects at high concentration, the accuracy of proton spectroscopic imaging could be improved. MR spectroscopic imaging is thus a unique tool to non-invasively measure the absolute temperature in vivo, but the low temporal resolution limit its applicability for real-time control of temperature except in case of long treatment duration.

Recently, reports have appeared about novel contrast agents with high temperature sensitivity. Two different mechanisms have been used to render the contrast temperaturesensitive. The first approach was to incorporate conventional contrast agents in temperaturesensitive liposomes [START_REF] Fossheim | Thermosensitive paramagnetic liposomes for temperature control during MR imaging-guided hyperthermia: in vitro feasibility studies[END_REF]. Such liposomes do not allow water to exchange between the liposome interior and bulk fluid. Above the transition temperature T c of the liposome, the membranes become "leaky" for water molecules, and the bulk fluid "feels" the relaxation effect of the contrast agent. The temperature sensitivity is therefore achieved via the exchange process. The second approach is based on a direct, temperature-driven transition between a diamagnetic state to a paramagnetic state using bi-stable molecular complexes [START_REF] Bartholet | Spin transition Molecular Materials: Intelligent Contrast Agents for Magnetic Resonance Thermometry[END_REF]. A qualitative summary of the principles and performance of the suggested methods is given in REF [START_REF] Quesson | Magnetic resonance temperature imaging for guidance of thermotherapy[END_REF].

2 Advanced temperature MRI

Rapid MR temperature imaging techniques

Optimal MR thermometry requires high temporal and spatial resolutions to precisely monitor temperature distribution within the targeted organ and to predict the effectiveness of thermal treatment. Ideally, temperature changes are monitored in real-time. In this paper, realtime is defined such that the update time is small compared to a significant temperature change or a temperature change that is potentially hazardous and thus allowing intervention.

For example, if the maximum temperature increase is less than 1°C per second, a "real-time" temperature mapping frequency of about 0.1 Hz is sufficient. If the deposited energy might lead to temperature changes of 10°C per second, a much higher temperature mapping frequency is required.

Temperature imaging of the abdomen and thorax is complicated due to complex organ displacement and deformation (mostly resulting from respiration and cardiac cycles), and the presence of flow in large vessels. In addition, magnetic susceptibility is usually spatially inhomogeneous (in particular near intestines, stomach, oesophagi, lungs,…), which may At present time, a number of studies based on the PRF technique, mainly performed at mid to high field (0.5T-1.5T) have demonstrated that precise temperature imaging can be obtained [START_REF] Weidensteiner | Stability of real-time MR temperature mapping in healthy and diseased human liver[END_REF][START_REF] Weidensteiner | Real-time MR temperature mapping of rabbit liver in vivo during thermal ablation[END_REF]. An example at the lower end of this range is shown in Figure 2.

The use of parallel imaging was successfully tested on human liver at 1.5T(37), as illustrated in Figure 3. Despite signal intensity loss inherent to parallel imaging, scan time reduction of a factor up to 2 could be achieved with a four-element coil, thus reducing artefacts linked to organ displacement. This acceleration of acquisition time can also be used to acquire a larger number of slices over the region of interest. In addition, such parallel imaging techniques could be very useful at high fields (typically 3T), taking advantage of the increase of signal intensity. In combination with segmented EPI acquisition techniques, parallel imaging can also be used to reduce the echo train length per repetition time at constant imaging time, with the advantage of reducing image distortions due to susceptibility effects.

The use of balanced fast field echoes sequences (SSFP, Balanced FFE, TrueFISP) was also investigated to derive temperature maps [START_REF] Paliwal | SSFP-based MR thermometry[END_REF][START_REF] Scheffler | Fast frequency mapping with balanced SSFP: theory and application to proton-resonance frequency shift thermometry[END_REF]. The advantage of this fast imaging technique is to offer high image quality with excellent anatomical detail. However, the PRF effect of the MR signal in such sequences is non-linear, and hard to separate from dynamic susceptibility effects. Its value for in vivo temperature mapping needs to be evaluated.

Despite MR thermometry was successfully tested on several organs at low (0.2T), mid (0.5T) and high field (≥1.5T), it remains necessary to improve volumetric coverage of the region of interest in mobile tissues, to precisely characterize local distribution of the temperature in the pathologic target and in surrounding healthy tissue. Novel techniques like slice tracking and real-time adaptive image techniques may become useful additions to the modern, PRF based, fast temperature imaging methods. Also, it is expected that the PRF methods will show improved performance at high field 3

Real-time image processing

As PRF-based temperature maps result from phase differences, this technique is prone to motion artefacts (patient motion, respiration, cardiac activity) [START_REF] De Senneville | Real-Time Artefact Corrections For Quantitative MR Temperature Mapping[END_REF]. When physical causes of the artifacts are reduced with the use of rapid imaging and image monitoring methods (see also section 3.3), image processing techniques can be used to further improve the quality of the thermometry. Processing of an image must be done fast to ensure real-time monitoring of temperature evolution (see also section 2.2). In practical terms, this implies that image processing must be done well within the update time between subsequent images.

Time scale of motion artifacts

Motion artefacts during MR thermometry sequence can be divided in two categories, intra-scan and inter-scan motion, based on the time scale of the motion as compared to the image acquisition time.

a) Intra-scan motion

This artifact is caused by the movement of an object during MR signal readout, yielding a low quality image, showing typical ghosting and blurring in the phase-encode direction. This artifact can be reduced by accelerating the imaging technique. However, a compromise on the acquisition parameters has to be made to obtain sufficient SNR while shortening acquisition duration. In the case of respiratory motion, typical acquisition duration of 400ms for one slice was found to be a reasonable value to obtain a sufficient SNR and to avoid intra-scan motion artifacts. On-line correction of this artifact is routinely performed during the image reconstruction process on state-of-the-art MRI instruments.

b) Inter-scan motion

Good image quality is obtained with moving organs but the region of interest is not correctly registered between successive images. This artifact can be generated by reflex or accidental patient motion and periodic motion like respiratory or cardiac activity.

Thermometry artifacts generated by motion of the susceptibility field

Although the spatial transformation of the reference phase image can be corrected, an unwanted phase shift cannot be fully suppressed. In fact, the local magnetic susceptibility is generally not fully uniform and the resulting phase variation ϕ ∆ can be described by the following equation ( 19) :

( ) E Z E T FT K K FT B T T B . . 3 1 . . . . . . 2 2 1 0 0         ∆         - + ∆ = ∆ - χ γ α γ ϕ , 0 ≠ K [9]
where ( )

Z Y X K K K K , , =
the position vector in the reciprocal space, and χ ∆ the modification of the susceptibility field. The other terms have been described above. The first term of the sum is due to the PRF-shift with temperature, whereas the second term arises when an object motion occurs. Both PRF-shift with temperature and movement of the patient lead to phase variation and it is not possible to separate the two contributions within a single phase image. Modeling of inhomogeneous field of susceptibility in-vivo is difficult and the sensitivity of PRF thermometry has limited its application to immobile organs and organs that can be easily restrained like muscle [START_REF] Hazle | Magnetic resonance imaging-guided focused ultrasound thermal therapy in experimental animal models: correlation of ablation volumes with pathology in rabbit muscle and VX2 tumors[END_REF], breast [START_REF] Hynynen | MR imaging-guided focused ultrasound surgery of fibroadenomas in the breast: a feasibility study[END_REF] and the prostate [START_REF] Peters | Magnetic resonance thermometry for predicting thermal damage: an application of interstitial laser coagulation in an in vivo canine prostate model[END_REF].

However, strategies allow correcting inter-scan motion related errors in PRF based MR thermometry, avoiding explicit modeling of the susceptibility field. Those techniques require first to estimate motion.

Motion reduction techniques

Motion restraining devices in the positioning of the patient are helpful. In addition, motion artefacts can be reduced significantly using rapid imaging techniques and navigator echoes. The simplest technique to reduce periodic motion artifacts is to synchronize the acquisition of the images to a stable period of the cycle (e.g. at the end of expiration for breathing). On line evolution of the respiration can be monitored with external methods using a pressure sensor positioned on the abdomen of the subject. Synchronization can be achieved by triggering the MR pulse sequence based on signals provided by the sensor ("respiratory gating" [START_REF] Morikawa | Feasibility of respiratory triggering for MR-guided microwave ablation of liver tumors under general anesthesia[END_REF]. Alternatively MR methods can be employed using so-called navigator echoes [START_REF] De Zwart | Online correction and visualization of motion during MRI-controlled hyperthermia[END_REF]). The major drawback of those methods is that the temporal resolution depends on the motion frequency and a limited set of images can be acquired during a respiratory cycle (~ 5 seconds for humans).

Physical methods for estimation of motion

Navigator echoes

Navigator echoes use an approach that consists of acquiring one or multiple lines of Fourier space including the center and applying an inverse Fourier transformation in order to obtain the profile of the object. This profile can be use to estimate motion [START_REF] De Zwart | Online correction and visualization of motion during MRI-controlled hyperthermia[END_REF]. However, the use of navigator echoes is restricted to rigid body motion and may not be optimal for complex organs displacements and/or deformations such as in the abdomen (liver, kidney) and the heart.

Ultrasonic echo

By sending an ultrasonic wave to a reflective target and measuring the time of flight for the return of this signal, the distance between the source and the target can be computed by multiplying this time by the speed of sound. This technique allows estimating the 3D translation of a point located in an unheated region. In general, the displacement of the heated region is hard to estimate in vivo because the coagulation caused by the heating modifies the echo. The technique is based on the hypothesis that the whole studied region undergoes the same displacement. Only the displacement of a point source of the target is obtained. This technique has been successfully used to target a mobile organ with an external heating source (for example with a focus ultrasound device) [START_REF] Pernot | 3D Real-time Motion Correction in High Intensity Focused Ultrasound Therapy[END_REF].

Active tracking coil

An active tracking coil can be used to obtain the 3D coordinates relative to the heat source location [START_REF] Suprijanto | Displacement Correction Scheme for MR-Guided Interstitial Laser Therapy[END_REF]. The active tracking coil is integrated into a power laser sheath that is introduced into the object. A custom receiver is used to detect location information from the active coil.

Like the preceding approach, only the displacement of a point source of the target is obtained.

Estimation of motion based on image registration techniques

This approach consists of using image processing techniques to detect organs displacements on the anatomic images. The objective is to relate the coordinate of each part of tissue in the image to register with the corresponding tissue in the reference image.

Contrary to the methods described above, complex organs displacements can be successfully estimated. To allow on-line monitoring of MR thermometry, this approach requires the use of a registration method that estimates motion on images under the following conditions

• no specification of the observed organ

• no user intervention

• real-time implementation.

In practice, under on-line conditions, it is often not practical to acquire isotropic 3D temperature maps on mobile organs, due to technical limitations of MR hardware and acquisition sequences. Thus, it is difficult to perform 3D image registration and one practical approach consists of estimating organ displacements in a 2D image generated from objects that are moving in 3D space [START_REF] De Senneville | Real-Time Artefact Corrections For Quantitative MR Temperature Mapping[END_REF]. Slices must be oriented in the principal axis of the organ displacement so that apparent movement is the best approximation of the real movement.

Image based registration can be divided in the following approaches :

Methods based on artificial foreign objects attached to the patient

Generally, such registration methods [START_REF] Lemieux | Effects of fiducial marker localization on stereotactic target coordinate calculation in CT slices and radiographs[END_REF] are difficult to be used for non-invasive MR therapy because :

• invasive markers cannot obviously be used [START_REF] Maurer | The accuracy of image-guided neurosurgery using implantable fiducial markers[END_REF][START_REF] Maurer | A review of medical image registration[END_REF].

• the nature of the registration transformation is often restricted to be rigid [START_REF] Maguire | Graphics applied to medical image registration[END_REF][START_REF] Leslie | External reference markers for the correction of head rotation in brain single-photon emission tomography[END_REF][START_REF] Simon | Accuracy validation in image-guided orthopaedic surgery[END_REF].

Methods based on patient generated image information only

Motion can be estimated on a set of landmarks, segmented structures, or computed directly from the image intensities (54):

-Methods based on landmarks may not be practical in our case because user interaction is required for the identification of the landmarks [START_REF] Evans | Anatomical-functional correlative analysis of the human brain using three dimensional imaging systems[END_REF] and the real-time control necessitates robust, and preferably automatic, methods.

-Segmentation based methods (with region extraction like deformable models [START_REF] Gueziec | Large deformable splines, crest lines and matching[END_REF], geometrical characteristic computation like high curvature pixels [START_REF] Davis | Contour-Based Motion Estimation[END_REF], or zeroing of the Laplacian [START_REF] Hildreth | The computation of the velocity field[END_REF], etcetera) could be used but, in practice, the segmentation techniques depend on the observed organ and the computation time is often not compatible with the real-time constraint. Moreover, the registration accuracy depends on the accuracy of the segmentation step.

-Methods based on image intensities without user prioritization or segmentation step.

Since some implementations of these algorithms are fast, and they do not require specifications with respect to the organ of interest nor any other user interaction, this approach is compatible with the conditions required to allow on-line monitoring of MR thermometry.

Two approaches may be followed:

- -Local transformation estimation using: a) block-Matching (the image to be registered is divided in blocks and the translation of each block is estimated in a research window [START_REF] Chen | A New Block-Matching Criterion for Motion Estimation and Its Implementation[END_REF]; it is supposed that the motion in the image is constant in the block; b) differential estimation methods of the optical flow (a velocity field is estimated assuming that intensity is conserved during displacement). A regularity constraint is also required, like 1) local regularity constraint that may optimize some local energy expression [START_REF] Lucas | An Iterative Image Registration Technique with an Application to Stereo Vision[END_REF] or 2) global regularity constraint which attempts to minimize a global energy function [START_REF] Schunck | Determining optical flow[END_REF]). This kind of method may be very useful for organs with complex motion and distortion like liver and heart.

Motion correction for on-line MR thermometry

Referenceless PRF Shift thermometry

This approach estimates the background phase from each acquired image phase. A polynomial fit to the background phase outside the heated region is performed using a weighted least-squares fit. Extrapolation of the polynomial to the heated region serves as the background phase estimate, which is then subtracted from the actual phase [START_REF] Rieke | Referenceless PRF shift thermometry[END_REF].

The major drawbacks of this approach are :

-The phase in the heated region is deduced from the phase outside the heated region. It is possible only if the susceptibility in the object is uniform. In addition, the phase φ has to be unwrapped. Phase discontinuities may be problematic.

-The method requires that phase perturbation induced by motion and by frequency drift during the intervention can be precisely expressed with a mathematical function (it depends on susceptibility in the observed organs and its neighbours but also of the actual motion).

-The method requires that enough pixels with sufficient SNR are located outside the heated region in order to allow a robust approximation of the mathematical function. In MR thermometry sequence, lipid suppression often generates images without sufficient signal in some locations. In addition, this method may be difficult to use when heated regions are close to interfaces between organs with different magnetic susceptibilities (lung/liver) or signal content (fat/water).

Correct spatial reference for thermometry

If a motion occurs between the instant 1 n t and n t , the relative temperature map after motion may be erroneous. A simple method consists of detecting motion vector to correct phase and temperature maps [START_REF] De Senneville | Correction Of Accidental Patient Motion For On-line MR Thermometry[END_REF]. Then the phase image acquired at n t ( n ϕ ) is taken as new phase image reference, and the temperature mapping at i t ( n i > ) is computed with :

( )k T T n i i . ' ϕ ϕ - + ∆ = ∆ [16]
where ' T ∆ is the n-1 temperature map after motion correction. Using this method, information about temperature changes between the instant 1 n t and n t is lost. Moreover, as mentioned previously [START_REF] De Zwart | Online correction and visualization of motion during MRI-controlled hyperthermia[END_REF], noise in the temperature map increases with the number of images corrupted by motion in a time series. In addition, uncertainty induced by the correction technique is propagated through all images obtained after each displacement. As a consequence, this technique can be applied when periodic motion has already been corrected with synchronization techniques.

Multi-baseline

A collection of multiple baseline reference images can be used to generate temperature maps. The selection of baseline image for temperature computation can be based on navigator echoes data, on the displacement of an active tracking coil [START_REF] Suprijanto | Displacement Correction Scheme for MR-Guided Interstitial Laser Therapy[END_REF], or on the displacement estimate on anatomical images with image registration algorithms [START_REF] De Senneville | Atlas-Based Motion Correction For On-Line MR Temperature Mapping[END_REF]. This approach can be used either to perform rapid acquisition (acquisition time per image faster than typical motion period) without synchronisation techniques in order to improve the temporal resolution [START_REF] Suprijanto | Displacement Correction Scheme for MR-Guided Interstitial Laser Therapy[END_REF][START_REF] De Senneville | Atlas-Based Motion Correction For On-Line MR Temperature Mapping[END_REF], or to improve the quality of the thermometry when periodical motion is reduced (69).

Image quality test

It is necessary to use a quality criterion to test the image processing efficiency. Indeed, intrascan motion and modification of the image plane position by through-plane motion can be detected rapidly. It is thus possible to detect those motion-corrupted images by using a realtime evaluation of a similarity criterion between the acquired anatomical images and a reference image [START_REF] Weidensteiner | Stability of real-time MR temperature mapping in healthy and diseased human liver[END_REF][START_REF] Suprijanto | Displacement Correction Scheme for MR-Guided Interstitial Laser Therapy[END_REF]. Images that do not match the quality criterion should not be included in the thermometry process.

Conclusion

The excellent linearity of the temperature sensitivity of the water PRF and its nearindependence with respect to tissue type, together with good temperature sensitivity, make PRF based temperature MRI the preferred choice for many applications, in particular at mid or high field (≥ 1 T), and probably even at very high field (≥ 3T). The PRF methods employ RF-spoiled gradient echo imaging methods in order to measure the phase change resulting (For complete details, see [START_REF] De Senneville | Atlas-Based Motion Correction For On-Line MR Temperature Mapping[END_REF]).

  induce important image distortions, especially in the case of rapid phase-sensitive, PRFbased, MR thermometry. A number of rapid imaging sequences, associated with specific reconstruction algorithms, have been developed for near real-time cardiac imaging and functional brain MRI (fMRI). These techniques for scan time reduction were based on the rapid acquisition of the full or a part of k-space (EPI imaging and associated segmented-EPI techniques), or SSFP (True FISP) imaging made possible with the hardware performance of modern scanners. More recently, the use of parallel imaging (eg SENSE[START_REF] Pruessmann | SENSE: sensitivity encoding for fast MRI[END_REF]) with multiple coils, UNFOLD[START_REF] Madore | Unaliasing by fourier-encoding the overlaps using the temporal dimension (UNFOLD), applied to cardiac imaging and fMRI[END_REF] have opened new possibilities to accelerate image acquisition. These new acquisition methods are very promising for the field or MR thermometry and they could be employed to partially overcome some difficulties of fast PRF thermometry, but with the cost of signal to noise reduction in most of the cases. As a consequence, an optimal compromise has to be made in the choice of the acquisition parameters for a given fast imaging sequence and a target organ, to avoid a dramatic increase of the uncertainty of phase estimate.PRF-based temperature MRI methods are based on gradient echoes with long TE (TE approximately equal to T2*). In conventional gradient-echo imaging sequences, the repetition time of the experiment is longer than the echo time. Therefore, a long TE for optimal sensitivity would require also a long TR, and temporal resolution is suboptimal. However, the use of echo shifting (leading to TE>TR in fast gradient echo sequences[START_REF] Moonen | A fast gradient-recalled MRI technique with increased sensitivity to dynamic susceptibility effects[END_REF]) avoids this compromise for MRI thermometry[START_REF] De Zwart | Fast lipid-suppressed MR temperature mapping with echo-shifted gradient-echo imaging and spectral-spatial excitation[END_REF][START_REF] De Zwart | Fast magnetic-resonance temperature imaging[END_REF]. So long as the length of the echo train is short compared with TE, a segmented EPI method can be used with multiple read-out periods per TR. When combined with echo-shifting (and maintaining RF spoiling) such a sequence is called a PRESTO sequence (PRinciples of Echo Shifting with a Train of Observation periods,[START_REF] Liu | A functional MRI technique combining principles of Echo-Shifting with a train of observations (PRESTO)[END_REF]). Therefore, PRF changes can be measured with very fast echo-shifted RF spoiled gradient echoes. Such methods show increased motion sensitivity, and may therefore be useful in immobile tissues. MRI sequences may be based on radial acquisition of the k-space instead of a standard rectilinear sampling, and thus offer the advantage of being less sensitive to motion. The principle is to acquire central points of k-space at each repetition time and to use these data as navigator echoes to correct each profile displacement.

  Global transformation estimation using: a) principal axes and moments based methods (rigid motions are estimated with the spatial distribution of their mass (59)); b) Fourier-based image registration algorithm (the properties of invariance in rotation and translation are used (60)); this method can register 2D images that are misaligned due to translation, rotation and scale (61)); c) based on a parameterized model of motion (62) (identification of the displacement is done by optimizing model parameters).

  from temperature-dependent changes in resonance frequency. When combined with modern parallel imaging techniques and other acceleration methods, the PRF based temperature maps can be acquired rapidly. A standard deviation of less than 1 o C for a temporal resolution below 1 s and a spatial resolution of about 2 mm is feasible for a single slice for immobile tissues. The rapid PRF methods generally allow only relative temperature measurements and fat suppression is necessary for most tissues. The methods require excellent registration to correct for displacements between scans since temperature maps are derived from relative changes between images. Navigator echoes and advanced image processing tools can be employed to detect and correct intra and inter-scan motion artifacts in real-time. Such tools will become indispensable for accurate guidance and control of local thermal therapies, in particular when treating mobile tissues.

Figure. 1 :

 1 Figure. 1: MR images showing the peak temperature increase after a 10-W sonication. Left:

Figure 2 .

 2 Figure 2. Isotherms following 6 minutes RF energy deposition in the paraspinal muscle of an

Figure 3 :

 3 Figure 3 : Transversal images in three liver tumor patients acquired with the respiratory-

Figure 4 :

 4 Figure 4 : example of motion correction in temperature mapping of the abdomen a a healthy

"Ligue contre le Cancer", France; the Conseil Regional d'Aquitaine, and