Baudouin Denis De Senneville

Mario Ries

Grégory Maclair

Chrit Moonen

MR-Guided Thermotherapy of Abdominal Organs Using a Robust PCA-Based Motion Descriptor

published or not. The documents may come   L'archive ouverte pluridisciplinaire

MR-guided thermotherapy of abdominal organs using a robust PCA-based motion descriptor Baudouin Denis de Senneville, Mario Ries, Grégory Maclair, and Chrit Moonen Abstract-Thermotherapies can now be guided in real-time using magnetic resonance imaging (MRI). This technique is rapidly gaining importance in interventional therapies for abdominal organs such as liver and kidney. An accurate online estimation and characterization of organ displacement is mandatory to prevent misregistration and correct for motion related thermometry artifacts. In addition, when the ablation is performed with an extracorporal heating device such as High Intensity Focused Ultrasound (HIFU), the continuous estimation of the organ displacement is the basis for the dynamic adjustment of the focal point position to track the targeted pathological tissue.

In this paper, we describe the use of an optimized PCA-based motion descriptor to characterize in real-time the complex organ deformation during the therapy. The PCA was used to detect, in a preparative learning step, spatio-temporal coherences in the motion of the targeted organ. During hyperthermia, incoherent motion patterns could be discarded, which enabled improvements in motion estimation robustness, the compensation of motion related errors in thermal maps, and the adjustment of the beam position.

The suggested method was evaluated for a moving phantom, and tested in-vivo in the kidney and the liver of twelve healthy volunteers under free breathing conditions. The ability to perform a MR-guided thermotherapy in-vivo during HIFU intervention was finally demonstrated on a porcine kidney.

Index Terms-Image registration, Motion analysis, Motion compensation, Magnetic resonance imaging

I. INTRODUCTION

M INIMALLY-invasive interventional procedures such as tissue ablation with radio-frequency (RF) and noninvasive ablation techniques such as high-intensity focused ultrasound (HIFU) show a high potential in oncology and cardiology as an alternative to classical surgery [START_REF] Cline | MR-guided focused ultrasound surgery[END_REF] [START_REF] Pernot | High power transcranial beam steering for ultrasonic brain therapy[END_REF]. The objective is a necrosis of the tissue caused by an elevated thermal dose due to sustained heating [START_REF] Sapareto | Thermal dose determination in cancer therapy[END_REF]. These methods are thus ideally paired with a non-invasive imaging modality which can provide both temperature measurements and target tracking information. The most promising candidate for this role is real-time MR-thermometry based on the water proton resonance frequency (PRF) [START_REF] Quesson | Magnetic resonance temperature imaging for guidance of thermotherapy[END_REF]. This technique provides continuous temperature mapping inside the human body as well as target tracking information by exploiting the complex nature of the MR-signal: whereas the signal magnitude provides anatomical information such as the position and the composition of the tissue, the phase ϕ is directly proportional Baudouin Denis de Senneville, Mario Ries and Chrit Moonen are with the Laboratory for Molecular and Functional Imaging: From Physiology to Therapy, FRE 3313 CNRS/University of Bordeaux 2, 33076 Bordeaux, France, (e-mail: {baudouin,gregory,ries,moonen}@imf.u-bordeaux2.fr).

Manuscript received ...; revised ... to the local magnetic field B and to the local proton resonance frequency (and thus the local temperature). The simplest way to obtain an estimate of the temperature changes at instant n (noted ∆T n ) is to evaluate phase shifts between dynamically acquired phase images and reference data sets as follows:

∆T n = (ϕ ref -ϕ n ) .k k = 1 γ.α. B.T E ( 1
)
where γ is the gyromagnetic ratio (≈ 42.58 MHz/Tesla), α the temperature coefficient (≈ 0.009 ppm/K) [START_REF] Peters | Ex vivo tissue-type independence in proton-resonance frequency shift mr thermometry[END_REF], and T E the echo time.

This new type of interventional procedure is very interesting for treatment of vital organs (kidney, liver, heart). However those organs are subject to physiological motion which leads to the following challenges:

1) Misregistration between phase images must be compensated, and temperature information must be mapped to a reference position of observed organs to obtain thermaldose measurements, requiring the voxel-by-voxel temporal integration of the temperature. 2) Reliable temperature measurement using PRF based MR-thermometry on abdominal organs is complicated by the fact that the target moves through an inhomogeneous and time-variant magnetic field, causing strong thermometry artifacts [START_REF] Poorter | The proton resonance frequency shift method compared with molecular diffusion for quantitative measurement of two dimensional time dependent temperature distribution in phantom[END_REF].

3) The estimation of the organ displacement is the basis for the dynamic adjustment of the HIFU beam to track the targeted pathological tissue, preventing an inefficient treatment and unwanted destruction of healthy tissue [START_REF] Denis De Senneville | Real time adaptive methods for treatment of mobile organs by MRI controlled high intensity focused ultrasound[END_REF]. To accomplish these tasks, the organ motion must be estimated and a number of techniques have been proposed in the past for various imaging modalities [START_REF] Maintz | A survey of medical image registration[END_REF] [START_REF] Mikic | Segmentation and tracking in echocardiographic sequences: Active contours guided by optical flow estimates[END_REF]. In the context of real-time guidance of the hyperthermia, navigator echoes [START_REF] De Zwart | On-line correction and visualization of motion during MRI-controlled hyperthermia[END_REF] or ultrasonic echoes have been originally proposed. However, only translational displacements can in those cases be estimated. Since modern MRI acquisition methods allow the rapid acquisition of large data volumes combined with an excellent tissue contrast and high spatial resolution, complex deformation can be estimated on a voxel-by-voxel basis using optical flow based approaches [START_REF] Roujol | Real-time MR-thermometry and dosimetry for interventional guidance on abdominal organs[END_REF].

Combined with the motion estimation process, several correction strategies have been proposed in the past to correct for thermometry artifacts, such as navigator echoes [START_REF] De Zwart | On-line correction and visualization of motion during MRI-controlled hyperthermia[END_REF] or referenceless phase corrections [START_REF] Denis De Senneville | Motion correction in MR thermometry of abdominal organs: A comparison of the referenceless vs the multi-baseline approach[END_REF]. An approach emerged recently which consists of analyzing phase perturbation with motion during a pre-treatment step performed prior to the intervention [START_REF] Denis De Senneville | Real time adaptive methods for treatment of mobile organs by MRI controlled high intensity focused ultrasound[END_REF]. For this purpose, a reference data-set of magnitude and phase images is recorded during the motion cycle of the organ. These images are used to evaluate a direct linear relation between the target motion and resulting phase variations.

Finally, since many applications in the field of MR-guided interventions also exploit the MR-image information for a direct retroactive control of the heating device [START_REF] Denis De Senneville | Real time adaptive methods for treatment of mobile organs by MRI controlled high intensity focused ultrasound[END_REF], both the update time and the latency of the obtained motion field are key points: for typical abdominal organ motion, it has been recently demonstrated that a minimal frame rate of 10 images/second is required and the latency (i.e the delay between the actual time of displacement and the availability of the motion information) must not exceed 100 milliseconds [START_REF] Ries | Real-time 3D target tracking in MRI guided focused ultrasound ablations in moving tissues[END_REF]. This latency is composed by the sum of the remaining acquisition time after echo-formation, the required data transport time, the image processing time, and the switching time of the HIFU-generator.

Black et al. proposed in 1997 a robust image-based motion descriptor (i.e. a small set of parameters representing the deformation) evaluated using a reduced parameterized flow models initially computed, during a preparative learning step, with a Principal Component Analysis (PCA) [START_REF] Black | Learning parameterized models of image motion[END_REF]. The technique was originally evaluated in the scope of human labial recognition. Recently, Mauldin et al. proposed a PCA-based displacement estimation in the context of ultrasound imaging [START_REF] Mauldin | Complex principal components for robust motion estimation[END_REF]. Previous work of the authors' group [START_REF] Maclair | PCA-based image registration : Application to on-line mr temperature monitoring of moving tissues[END_REF] [START_REF]PCA-based magnetic field modeling : Application for on-line mr temperature monitoring[END_REF] explored the possibility to apply the PCA-based motion descriptor for real-time MR-guided thermotherapy of mobile organs. The PCA-based motion descriptor was used to detect, in a preparative learning step, spatio-temporal coherences of the periodic organ motion. During hyperthermia, incoherent motion patterns could be discarded, which allowed the following improvements:

1) Optical-flow based algorithms rely on the assumption of conservation of local intensity along the trajectory which can be violated during thermotherapy because rapid MRimaging is in general associated with low Signal-To-Noise ratio (SNR). In addition, since the tissue is heated, several MR relevant tissue properties such as T 1 , T 2 and T * 2 relaxation times are subject to change during imaging [START_REF] Graham | Time and temperature dependence of MR parameters during thermal coagulation of ex vivo rabbit muscle[END_REF]. This leads to local intensity variations, which in turn can be misinterpreted by optical-flow based algorithms as "motion". The PCA-based motion descriptor was used to provide a flow field that was consistent with the learned model and robust under the assumption of global brightness constancy but allowed local intensity variations.

2) The PCA-based motion descriptor was used to model the magnetic field variation with the target displacement. This improved the correction of motion related errors in temperature maps.

In the current paper, this PCA-based method was modified as follows: We first briefly introduce a method, based on the PCA-based representation, which is designed to separate physiological displacement from noise contributions. This step was found to be a necessary prerequisite for the practical application of the method. The potential applications enabled by the technique are then described and tested with a phantom experiment, and an in-vivo thermometry stability study in the kidney and the liver of twelve healthy volunteers under free breathing conditions. The ability to perform a MR-guided thermotherapy in-vivo during a real intervention was finally demonstrated on a porcine kidney during a HIFU heating experiment.

II. METHOD DESCRIPTION

During the therapy, MRI is used for two tasks: MR thermometry and target tracking. Since the basis of the PCAbased approach is described in detail in [START_REF] Maclair | PCA-based image registration : Application to on-line mr temperature monitoring of moving tissues[END_REF] [START_REF]PCA-based magnetic field modeling : Application for on-line mr temperature monitoring[END_REF], only a short summary is given here. Details are given for the method designed to separate physiological displacement from noise contributions, based on the PCA-based representation.

We evaluated the method in 2D since it is in practice possible for respiratory induced motion, to choose the imaging plane direction parallel to the principal axis of the organ displacement.

A. Preparative learning step

Motion patterns were learned during a preparative learning step performed before hyperthermia, based on a training set of N images (magnitude and phase). The motion cycle was sampled with a sufficient density in order to avoid discretization errors. With a typical imaging frame rate of 10 images per second and a respiration frequency of 3-5 seconds this pretreatment step can be completed in a relatively short duration of 15-20 seconds (in this study, we used N = 200 images).

1) Training set of organ displacement: A model of image motion was obtained from a training set of N flow fields. For that purpose, an affine transformation was first estimated using a differential Gauss-Newton approach [START_REF] Friston | Spatial registration and normalisation of images[END_REF]. The obtained motion field was used for preconditioning of a more complex deformation estimation: The implemented approach was originally proposed by Horn&Schunck [START_REF] Horn | Determining optical flow[END_REF] for optical image stabilization and computed a displacement field on a pixel-by-pixel basis. It includes a modification proposed by Cornelius&Kanade [START_REF] Cornelius | Adapting optical flow to measure object motion in reflectance and x-ray image sequences[END_REF], which relaxes the condition of intensity conservation from a global to a regional conservation. Finally, in order to optimize the computation time and to stabilize the convergence of the algorithm, a multi-resolution scheme was used [START_REF] Pratikakis | Robust multiscale deformable registration of 3D ultrasound images[END_REF] which iterated the registration algorithm from a four-fold down-sampled image step-by-step to the full image resolution.

2) Learning parameterized flow models: At this point, a model of image motion was obtained, which was a set of N optical flow fields. The set was used to build a parameterized flow model. PCA was used to find an orthonormal basis that spans an N -dimensional vector space and depict the underlying characteristic patterns of the motion cycle.

3) Estimation of the motion descriptor: The periodic motion cycle was characterized with a small set of parameters as follows: The approximated spatial transformation T t was, for each pixel (x, y) calculated as a linear combination of the previously computed basis vectors B i :

T t (x, y) = M -1 i=0 D t i B i (x, y) (2)
where D t i (0 ≤ i < M ) was the motion descriptor (i.e. a set of parameters representing organ displacement), M (M ≤ N ) the size of the motion descriptor, and 0 ≤ t < N was the index of each image of the training set. The objective was thus to find the coefficients D t i that produced a flow field minimizing the matching error:

LS = M ref -T -1 t (M t ) 2 (3)
where M ref was the reference magnitude image and M t the actual image. This minimization was realized using a Marquardt-Levenberg least square (LS) solver [START_REF] Marquardt | An algorithm for least squares estimation on nonlinear parameters[END_REF]. To stabilize the convergence of the algorithm, a spatial Gaussian filter (kernel 3×3) was applied on the new incoming image M t , and the motion descriptor obtained from the previous acquisition was used as a starting point for the new motion estimation.

4) Separation of physiological displacements from noise contributions:

We used the PCA-based representation of the motion to detect spatio-temporal coherences of the periodic organ motion, in order to separate physiological displacements from noise contributions. For that purpose, only the M eigenvectors B i associated with the m largest eigenvalues λ i (0 ≤ i < M < N ) were conserved as follows. The objective was thus to determine the threshold value for M , which separates eigenvectors representing physiological motions (i ≤ M ) and from eigenvectors coding for noise contribution (i > M ). To accomplish this task, possible values for M were iteratively enumerated from 2 to N (we note m the iteration number), and, for each m, the following test was performed:

1) The PCA-based motion descriptors D t i were computed for all acquisitions of the preparative learning step.

2) A temporal variation of the m th coefficient of the motion descriptor (i.e. D t m-1 ) was analyzed over the preparative learning step, as follows. Since the respiratory and the cardiac activities are periodic, we analyzed D t m-1 , ∀t, 0 ≤ t < N in the Fourier domain. Typical periods of the respiratory and cardiac activities are, in the general case, in the range of 3 -6 seconds and 0.5 -2 seconds, respectively. A threshold of 4 Hz was used to separate physiological motions from the noise contribution, and the maximal spectral magnitude was compared for lower and higher frequencies. The corresponding (m -1) th eigenvector was considered to be noise corrupted if the spectral magnitude for frequencies above 4 Hz exceeds the value of 20 % of the main peak below 4 Hz. M was then set to (m -1), and the iterative enumeration was stopped.

5) Modelling motion induced phase changes:

The N images of the preparative learning step were used to evaluate the linear relation between organ motion and phase variation. For that purpose, the obtained motion field was used to remap all pixels of the phase images to their reference position. The overall magnetic field variations were approximated as the sum of linear phase changes of each principal motion component on a pixel-by-pixel basis:

T -1 t (ϕ t (x, y)) = M -1 i=0 D t i P i (x, y) + P M (x, y), 0 ≤ t < N (4
) where P i (0 ≤ i < M ) denotes the parameterized magnetic field model, with P M representing the initial phase distribution, and D t i were the coefficients calculated from Eq. ( 2) and (3). From the set of N equations with m unknowns defined by ( 4), an overestimated system could be obtained if m N . P i (x, y) are individually computed for each pixel (x, y) using a Singular Value Decomposition (SVD).

B. Hyperthermia procedure

During hyperthermia, each time a new image was acquired, the PCA-based motion descriptor was estimated and used to compute (with Eq. ( 3)) a flow field that was consistent with the learned model. The PCA-based motion descriptor was then used to correct for motion related errors on temperature maps as follows: During hyperthermia, the parameterized magnetic field model P i allows the reconstruction of magnetic field distribution corresponding to the current position of an organ. A synthetized reference non-heated phase image was calculated using the motion descriptor D t i representing the actual organ displacement as follows:

ϕ ref (x, y) = M -1 i=0 D t i P i (x, y) + P M (x, y) (5)
In this study, we set an upper bound of 6 for M , in order to keep the computational time for the motion estimation process of one image below 25 ms with the used experimental setup.

The measured registered phase image was subtracted to the reference phase map ϕ ref to compute a motion corrected temperature map (with Eq. ( 1)). Since the 2π periodicity of the image phase could lead to aliasing artifacts in the temperature maps, a temporal phase unwrapping on a pixel-by-pixel basis was applied. This process was valid under the condition that the temperature variation between two successive acquisitions did not create a phase variation greater than π.

C. Experimental setup

Dynamic MR temperature imaging was performed on a Philips Achieva 1.5 T (Philips Healthcare, Best, The Netherlands). MR images were processed using a dual processor (INTEL 3.1 GHz Penryn, four cores) workstation with 8 GB of RAM and dual 1 GB/s network interface cards.

Three experiments were performed to assess the improvement obtained with the PCA-based motion descriptor:

1) The method was applied on an ex-vivo heated phantom prone to a periodical translational displacement. The estimated motion was compared to reference positions provided by an external sensor.

2) The potential of the proposed method to modelize the magnetic field perturbation and to improve the correction for motion related errors on thermal maps was evaluated on the abdomen of 12 healthy human volunteers under free-breathing conditions. 3) A real-time tracked HIFU heating experiment was performed on a porcine kidney. Since the kidney displacement was rigid, the target motion could be succesfully estimated on-line by a translation restrained on a ROI manually set on the kidney at the begining of the hyperthermia (a differential Gauss-Newton approach was used [START_REF] Friston | Spatial registration and normalisation of images[END_REF]). The estimated displacements were subsequently used as reference to evaluate off-line the proposed PCAbased method. It is important to see those proposed experiments in the appropriate context: Our PCA-based method presents a real interest only for elastic organ deformation such as in the liver, and was not a necessary prerequisite for experiments #1 and #3. However, those two experiments were designed to quantify, using reference positions, the accuracy and the precision of the proposed PCA-based method under thermoablation conditions.

1) Ex-vivo heating study: 3000 images were dynamically acquired during 3 minutes using a dual-shot, gradient recalled echo-planar acquisition sequence with the following parameters: repetition time (T R)=30 ms, T E=15 ms, bandwidth in readout direction=1777 Hz,flip angle=20 • , field of view (F OV )=256×104 mm 2 , slice thickness=5 mm, matrix=128× 58. A porcine muscle was positioned on a motorized plateform, which generated a periodic translational displacement (amplitude=10 pixels and frequency=0.5 Hz). RF heating was performed 10 seconds after the end of the preparative learning step, using a clinical MR-compatible bipolar RF device (Radionics, Burlington, MA) with 8 W of RF-power during 75 seconds. For an independent assessment of the object displacement, an additional navigator echo, which provided a one dimensional displacement information (0.5 mm precision), was positioned parallel to the displacement on the apex of the muscle in order to get the reference displacement. We compared temperature maps obtained when phase images of Eq. ( 5) were registered with the navigator echo (which serves as reference results), with the Cornelius&Kanade algorithm (the algorithm is described in details in section II.A.1) and with the PCA-based motion descriptor.

2) On-line in-vivo study: Dynamic MRI was performed under free-breathing conditions on the abdomen of 12 healthy human volunteers under real-time conditions. An imaging frame rate of 10 images/s was maintained for five minutes of MR-imaging while MR-Thermometry was performed in real-time. The MR sequence was a single-shot gradient recalled echo-planar which employed the following parameters: 3000 dynamic coronal images, one slice, T R=100 ms, T E=26 ms, bandwidth in readout direction=2085 Hz, flip angle=35 • , F OV =256 × 168 mm 2 , slice thickness=6 mm, matrix=128 × 84, using a four element phased array body coil. No hyperthermia was performed and the statistical evaluation of the temperature stability was performed on the kidney and the liver of each volunteer individually by averaging the temporal temperature standard deviation over a manually defined ROI. The thermometry was compared with a synthetized reference non-heated phase expressed in Eq. ( 4) and ( 5), as the linear combination of the six free parameters of an affine transformation [START_REF] Roujol | Real-time MR-thermometry and dosimetry for interventional guidance on abdominal organs[END_REF] (a differntial Gauss-newton approach was used for that purpose).

3) In-vivo heating study on a porcine kidney: MRI guided HIFU was performed in vivo in the kidney of a pig under general anesthesia and artificial breething. An imaging frame rate of 10 images/s was maintained for 2.5 minutes using a single-shot gradient recalled echo-planar sequence. The MR sequence employed the following parameters: 1500 dynamic sagittal images, one slice, T R=100 ms, T E=41 ms, bandwidth in readout direction=2085 Hz, flip angle=35 • , F OV =320×140 mm 2 , slice thickness=6 mm, matrix=128 × 56, using the integrated three elements phased array coil of the HIFU system. A MR compatible Philips HIFU ablation system (Sonalleve, Philips Healthcare, Helskinki, Finland) composed of a table top containing a 256 elements HIFU transducer, integrated in the 1.5 T Achieva-Intera MRI was used to perform a temperature elevation. The transducer radius and aperture were 120 mm and 126 mm, respectively, creating an ellipsoid focal point (1×1×7 mm 3 ). The animals were placed in right lateral decubitus position so that the kidney was accessible through an unobstructed beam-path directly below the rib-cage. Motion compensated MR-thermometry and tracking was performed for a duration of 1 minutes on the kidney with power of 100 W. A temporal filtering based on an infinite impulse-response filter was applied to the temperature maps to increase the SNR (more details can be found in [START_REF] Roujol | Real-time MR-thermometry and dosimetry for interventional guidance on abdominal organs[END_REF]).

III. RESULTS

A. Ex-vivo heating experiment

Fig. 1 shows thermal maps obtained after 50 seconds of RF heating obtained using the reference navigator echo based registration (1a), the Cornelius&Kanade algorithm (1b) (note the deformation of the heated area), and the PCA-based motion descriptor are displayed (1c). A detail of the anatomical image in the heated area before hyperthermia and an image obtained after 50 seconds of RF heating are given in Fig. 1a. The position of each electrode of the bipolar RF device is indicated with the two red arrows. A strong signal decrease can be observed during the hyperthermia process in the heated area (see green arrow) due to a large temperature change. It can be observed that the PCA-based motion descriptor provided temperature maps which were similar to the reference results.

An analysis of the estimated motion field accuracy over the whole duration of the experiment is presented in Fig. 2a. For that purpose, the euclidean error between the estimated displacement and the navigator echo information was computed in a pixel located in the heated area (reported by the green arrow in Fig. 1a). It can be observed in the first 10 seconds (i.e. before the energy deposition) that the precision was comparable with the PCA-based registration. However, the Cornelius&Kanade approach provided poor results during hyperthermia procedure starts (up to 6 millimeters of error), while the PCA-based method remained stable and accurate over the time (less than 2 millimeters of error). (2b), and with the PCA-based registration (2c). Compared to the reference navigator echo results (red line), the temperature measured using the Cornelius&Kanade registration was overestimated during the heating process. Contrary to the implemented Cornelius&Kanade approach, the temperature evolution obtained with the PCA-based motion descriptor was close to the reference results obtained with the navigator echo.

B. On-line in-vivo study

The volunteers studies depicted various conditions: Over the twelve volunteers, the SNR was evaluated as 10 ± 2.5 (min=7, max=14) in the kidney and 7 ± 3 (min=4, max=14) in the liver. The peak-to-peak motion amplitude obtained from the landmark points was 10 mm ± 4.5 (minimum=4, maximum=18) in the kidney and 11 mm ± 4.5 (minimum=6, maximum=18) in the liver. Due to variations of the respiratory volume, the peak-to-peak amplitude increased between the preparative learning step and the interventional procedure in patients #8, #9 and #12 (where an increase of 50 to 65 % was measured). This value exceeded 150 % in patient #7. Less than 10 % of increase was measured in patients #2, #3, #5 and #10.

Fig. 3 details the improvement of the temperature precision when the PCA-based motion descriptor was used to characterize the organ deformation in Eq. ( 4) and ( 5), as compared to an affine motion model. In average over all volunteers, the temperature stability was improved from an initial value of over 2.3 • C to 2 • C (kidney) and 3.3 • C to 2.7 • C (liver) using the PCA-based motion descriptor. None of the data were rejected due to misregistration or faulty phase correction due to the absence of spontaneous motion.

Fig. 4 and 5 detail results obtained in the volunteer for which the most significant motion amplitude variation was detected between the preparative learning step and the hyperthermia procedure (i.e. volunteer #7). Fig. 4 shows the vertical displacement estimated during the preparative learning step (4a) and during hyperthermia (4b) using the PCA-based registration (solid line). Manually tracked results provided by a staff scientist with a precision of a pixel (red dots) in a landmark located in the liver are also reported for comparison (red arrow in Fig. 5c). Due to deep variations of the breething, an increase of the peak-to-peak amplitude by 150 % can be observed between the preparative learning step and the interventional procedure in this case. Even in this case, the proposed correction ensured 2 • C of temperature stability in 70 % of all pixels of both the kidney and the liver, as shown in Fig. 5b. The first and the fourth eigenvectors are displayed in Fig. 5c and5d, respectively (for each pixel the amplitude of the displacement vector was computed. Each map was normalized between 0 and 1 for an easier vizualisation). Each of these eigenvectors represented a physiological component of the organ deformation: while the first component modelized the respiratory motion (motion field vectors were oriented in the head-foot direction), the fourth component represent artereal pulsations caused by the cardiac cycle. The estimation of a complex deformation was mandatory in the liver where the motion amplitude varied with a factor two (see Fig. 5c). C. In-vivo heating study on a porcine kidney Fig. 6 shows an overlay of the temperature map after 60 s of sonication with 100 W accoustic power without and with adjustment of the HIFU focal point position with respect to the target position. A low SNR value of 6 was obtained due to the difficulty to perform an accurate MR receiver coil placement. The tracking algorithm (which used a translational motion model) remained locked-on the observed motion pattern over the entire duration of the sonication experiment without apparent failure. The temperature evolution The estimated displacement was used off-line as reference to quantify the quality of the proposed PCA-registration. A signal decrease, correlated with the temperature rise could be observed in the heated area (see the time curve over the whole experiment of the estimated vertical displacement in a voxel located in the heated area of Fig. 7a). It can be observed in Fig. 7a and 7b that only the first eigenvector visually matched properly the translational kidney displacement. While the first component of the motion descriptor visually described properly the motion pattern induced by the assisted respiration (see Fig. 7d), the second component was very noisy (see Fig. 7f). The period of the respiratory cycle was ≈ 5 seconds and the main peak located at 0.2 Hz on the spectral magnitude of the first component (see 7e). The ratio between the maximal spectral magnitude for frequencies lower and higher than 4 Hz was 0.03 for the first component, and reached 0.69 for the second (see Fig. 7e and7g). The optimized size of the motion descriptor was thus of 1 component. Similarly to the ex-vivo experiment, the Cornelius&Kanade approach shows a shift of the estimated displacement correlated with the magnitude signal decrease (see Fig. 7.a and 7.h). Using the reference translation displacement estimated on-line, we evaluated that this shift reached a maximal value of 1.5 millimeters at the end of the heating process. The PCA-based approach remained stable over the time since less than 0.1 millimeters of error was estimated.

IV. DISCUSSION

The introduction of a PCA-based motion descriptor allowed several improvements compared to the previously established methods described in [START_REF] Roujol | Real-time MR-thermometry and dosimetry for interventional guidance on abdominal organs[END_REF] [START_REF] Ries | Real-time 3D target tracking in MRI guided focused ultrasound ablations in moving tissues[END_REF]. First, the method allowed to maintain only physiological components of the estimated motion in the registration process, thus allowing:

1) To reduce the noise on the estimated displacement (as shown in the ex-vivo study (Fig. 2.a) and the in-vivo heating experiment (Fig. 7.i)). 2) To discard during the hyperthermia procedure, possible local intensity perturbations induced by the heating: While purely optical-flow-based realignment may lead to temperature map computation errors for the case of local or global intensity changes, PCA-based realignment can give accurately registered temperature maps for the case of periodical motion, since it relies on a global fit of the principal components. Since the proposed method employs optical flow only during the preparative step and relies on a global fit of the image content during the intervention, estimated transformation T t was thus less susceptible to those local intensity variations. It is important to note that the proposed motion estimation was not constrained to estimate positions present in the collection. The proposed approach assumes that the deformation must be properly described by a linear combination of the flow basis P i . In the tested data, the proposed motion model was shown to be efficient to estimate organ motion amplitudes higher than the ones observed in the learning phase (see Fig. 4).

In addition, the method provided a characterization of the organ displacement which allowed, in turn, an accurate modelling of the magnetic field perturbation. As a result, the PCA-based motion descriptor improved the precision of the temperature computation, compared to an affine motion model. The proposed approach also allowed to extrapolate motion amplitudes higher than the ones observed in the preparative learning step and could still provide an estimate of the reference phase in those cases. It is interesting to note that few temperature precision improvements were obtained on patients depicting small motion amplitude variation between the preparative learning step and the hyperthermia procedure (see volunteers #2, #3, #5 and #10 in Fig. 3). Although an affine motion model was found to perform well in many patients, the PCA-based approach should be preferred to guarantee good treatment in such outlier cases. The method provided thermal maps registered to a reference position, which is a necessary prerequisite for the computation of the thermal dose [START_REF] Sapareto | Thermal dose determination in cancer therapy[END_REF], and thus for the on-line assessment of the progress of the interventional procedure. It must be emphasized that, in the current study, the magnetic field perturbation was estimated assuming a simple linear magnetic field variation with organ displacement. Although this assumption holds in general for small displacements, the precision of this simple model must be carefully evaluated in regions displaying large motion amplitudes or large susceptibility variations, such as in the vicinity of the digestive tubes (see arrow (i) in Fig. 3c), of the quadratus lumborum muscle (see arrow (ii) in Fig. 3c), of the vertebral column (see arrow (iii) in Fig. 3c), and in the upper part of the liver (see arrow (iv) in Fig. 3c). This can in practice be achieved by mapping the fitting error from Eq. ( 4). This also allows to discard regions where low signal levels in conjunction with large susceptibility variations prevent to adjust successfully the linear model to the phase data. The determination of the size M of the subset of principal components was a key point. In all presented experiments, the energy of the first eigenvalue was higher than 99 % of the total energy. This can be attributed to the fact that motion of abdominal organs is dominated by linear shifts due to the displacement of the diaphragm during the respiratory cycle. The determination of M could thus not be based on a pure statistical analysis. The proposed criterion based on the Fourier representation of each component of the motion descriptor allowed to discriminate physiological and noise contributions: An optimal size of 1 was determined for the in-vivo porcine kidney experiment which matched the kidney displacement with the used imaging plane orientation. A minimal size of 5 was required in the human study, since the liver depicts a complex deformation together with arterial pulsations. Althought a size higher than 6 would have been set in 6 cases with the proposed criterion (note that this provided improved thermometry results in Fig. 3), we limited this value to 6 for computational time reasons.

Ries et al. have shown that, for typical abdominal organ motion, a temperature elevation comparable to static control experiments can be achieved for update frequencies of more than 10 Hz and latencies of 100 ms [START_REF] Ries | Real-time 3D target tracking in MRI guided focused ultrasound ablations in moving tissues[END_REF]. On our test platform, the latency was composed as the sum of: the remaining acquisition time after echo-formation (≈50 ms), the required data transport time (≈10 ms), the image processing time, and the switching time of the HIFU-generator (≈10 ms). Thus, computation time for the whole image processing must typically be below 30 ms to be compatible with the real-time constraint. The PCA-based image realignment has computational advantages: the reduction of complex periodic motion patterns to the most significant principal components minimizes the degrees of freedom for the registration without a priori assumptions or simplifications of the form of the motion. On our implementation, the computation time required for the processing of one image was below 30 ms. Depending on the displacement amplitude of the target, between 8 and 25 ms were required for the motion estimation process. Less than 1 ms was necessary for the calculation of the motion corrected temperature map. This indicates that the proposed method is suitable for real-time interventional MR-guidance, as shown in the in-vivo heating study.

Although the proposed method improves previously suggested approaches for MR-guided thermotherapy of mobile organs, several limitations remain. The first limitation is the additional scan required to learn the motion pattern before hyperthermia. However, in practice this limitation is less severe, since the proposed correction approach is paired with high frame rate imaging and thus a re-calibration can be completed in the relatively short time of two to five respiratory cycles. The second problem arises when new positions are observed during hyperthermia, for which the learned motion pattern is inadequate. In this case, the proposed approach must be paired with a correction adapted for accidental motion, as detailed in [START_REF] Denis De Senneville | Motion correction in MR thermometry of abdominal organs: A comparison of the referenceless vs the multi-baseline approach[END_REF]. Finally, the effect of through plane motion is a severely limiting factor. In the scope of the presented study, this has been avoided by aligning the slice parallel to the main direction of motion. Since real-time acquired MR-images have generally an anisotropic image resolution (roughly 2-3 times higher in the slice than perpendicular to the slice) this did in none of the volunteers lead to severe problems due to residual through-slice motion. However, this is an inconvenient limitation of the degree of freedom for slice positioning (which for diagnostic reasons might be unfavorable) and complicates the slice placement in the planning of the intervention. Furthermore, it can not be excluded that patients with liver pathologies display more complex 3D motion patterns than the examined cases. In such cases the only true remedy would be to acquire extended 3D data volumes and apply the proposed methods in 3D. However, although it is well established that MR-imaging can provide motion estimates with a high spatial resolution, it is difficult to acquire on-line 3D isotropic images because technical limitations lead to a compromize between spatial and temporal resolution and low SNR associated with fast 3D acquisition sequences. Alternative approaches, based on the use of additional informations such as navigator echoes, in combination with adaptive slice tracking to address 3D trajectories, may be used [START_REF] Ries | Real-time 3D target tracking in MRI guided focused ultrasound ablations in moving tissues[END_REF].

V. CONCLUSION

The proposed approach addresses both motion compensated MR thermometry and target tracking by applying high frame rate MRI coupled with a real-time motion estimation and characterization obtained from all incoming images. The PCAbased motion descriptor introduced in [START_REF] Maclair | PCA-based image registration : Application to on-line mr temperature monitoring of moving tissues[END_REF] and [START_REF]PCA-based magnetic field modeling : Application for on-line mr temperature monitoring[END_REF] presented several advantages compared to previously suggested approaches: the image-based estimated motion was robust with respect to local and global intensity changes, and artifact-free temperature maps were obtained in all examined cases. The method allowed to achieve a sub-second temporal resolution with very short image latencies over sustained imaging periods of several minutes. During the intervention both, the target location and the target temperature were continuously available with a high temporal resolution and precision.

This renders the method well suitable for the MR-guidance of a heating intervention on abdominal organs in vivo under free-breathing over sustained periods of several minutes and presents therefore a step towards clinical non-invasive HIFU therapies of kidney and liver tumors.

Fig. 1 :

 1 Fig.1shows thermal maps obtained after 50 seconds of RF heating obtained using the reference navigator echo based registration (1a), the Cornelius&Kanade algorithm (1b) (note the deformation of the heated area), and the PCA-based motion descriptor are displayed (1c). A detail of the anatomical image in the heated area before hyperthermia and an image obtained after 50 seconds of RF heating are given in Fig.1a. The position of each electrode of the bipolar RF device is indicated with the two red arrows. A strong signal decrease can be observed during the hyperthermia process in the heated area (see green arrow) due to a large temperature change. It can be observed that the PCA-based motion descriptor provided temperature maps which were similar to the reference results.An analysis of the estimated motion field accuracy over the whole duration of the experiment is presented in Fig.2a. For that purpose, the euclidean error between the estimated displacement and the navigator echo information was computed in a pixel located in the heated area (reported by the green arrow in Fig.1a). It can be observed in the first 10 seconds (i.e. before the energy deposition) that the precision was comparable with the PCA-based registration. However, the Cornelius&Kanade approach provided poor results during hyperthermia procedure starts (up to 6 millimeters of error), while the PCA-based method remained stable and accurate over the time (less than 2 millimeters of error). Fig.2also shows the temporal evolution of the temperature in the same pixel with the Cornelius&Kanade registration method

Fig. 2 :

 2 Fig. 2: (a): Temporal analysis of the registration accuracy in a pixel located in the heated area (green arrow on Fig. 1a) obtained with the Cornelius&Kanade algorithm (black line) and the PCA-based motion descriptor (green line). (b): Temporal evolution of the temperature obtained with phase images registered with the Cornelius&Kanade algorithm. (c): Temporal evolution of the temperature obtained with phase images registered with the PCA-based motion descriptor. The reference temperature obtained with navigator echo based displacements is displayed with the red line in (b) and (c).

Fig. 3 :Fig. 4 :

 34 Fig. 3: Histogram of the temperature precision obtained over the twelve volunteers in the kidney (a) and the liver (b), when the temperature maps were computed using an affine motion model (dash line) and using the PCA-based motion descriptor (solid line).

Fig. 5 :

 5 Fig. 5: Results obtained in the abdomen of a healthy volunteer (volunteer #7 in Fig. 3): (a) and (b) show the temperature standard deviation map obtained with an affine motion model and the PCA-based motion descriptor, respectively. (c) and (d) display two PCA components (B 0 and B 3 ): For each pixel the amplitude of the displacement vector was computed, and each map was normalized between 0 and 1 for vizualisation. Results are presented in superimposition with the anatomical image depicting the kidney (arrow (1)) and the liver (arrow (2)).

Fig. 6 :

 6 Fig. 6: Thermometry results obtained for the in-vivo HIFU heating experiment. T2*-weighted magnitude images of the kidney of the untracked (a) and the tracked HIFU-ablation experiment (b) are overlayed with the temperature maps obtained after 60 s of sonication. The temporal evolution of the temperature in a voxel located in the heated area is reported in (c) for the untracked (black dashed line) and the tracked (red line) HIFU-ablation experiment. The untracked experiment reached a lower final temperature than the experiment performed with full motion compensation (tracked: 15 • C; untracked: 10 • C).

Fig. 7 :

 7 Fig. 7: Motion estimation results otained for the in-vivo HIFU heating experiment. The temporal evolution of the magnitude signal in the voxel with the maximal temperature rise is shown in (a) (note the signal intensity loss during the heating process). The two first PCA components (B 0 and B 1 ) are reported in (b) and (c), respectively. The temporal evolution of the two first coefficients of the PCA-base motion descriptor (i.e. D 0 and D 1 ) and their corresponding representation in the spectral domain are displayed in (d), (e), (f) and (g). The red dashed line in (e) and (g) shows the value of 20 % of the main peak below 4 Hz. The temporal evolution of the estimated displacement in the head-foot direction with the Cornelius&Kanade algorithm and the PCA-based approach are displayed in (h) and (i), respectively.

ACKNOWLEDGMENT

The authors would like to thank the "Ligue Nationale Contre le Cancer", the "Conseil Régional d'Aquitaine", "Diagnostic Molecular Imaging" EC-FP6-project LSHB-CT-2005-512146, the "Agence Nationale de la Recherche" (project "MRgHIFU-ALKT"), Foundation InNaBioSanté (project ULTRAFITT) and Philips Healthcare.


Sébastien Roujol

Mario Ries

Chrit Moonen

Baudouin Denis De Senneville

Automatic Non-Rigid Calibration of Image Registration for Real Time MR-guided HIFU ablations of mobile organs

come   L'archive ouverte pluridisciplinaire

Automatic Non-Rigid Calibration of Image Registration for Real Time MR-guided HIFU ablations of mobile organs Sébastien Roujol, Mario Ries, Chrit Moonen, and Baudouin Denis de Senneville

Abstract-Real time magnetic resonance imaging (MRI) is rapidly gaining importance in interventional therapies. An accurate motion estimation is required for mobile targets, and can be conveniently addressed using an image registration algorithm. Since the adaptation of the control parameters of the algorithm depends on the application (targeted organ, location of the tumor, slice orientation, ...), typically an individual calibration is required. However, the assessment of the estimated motion accuracy is difficult since the real target motion is unknown.

In this paper, existing criteria based only on anatomical image similarity are demonstrated to be inadequate. A new criterion is introduced, which is based on the local magnetic field distribution. The proposed criterion was used to assess, during a preparative calibration step, the optimal configuration of an image registration algorithm derived from the Horn and Schunck method.

The accuracy of the proposed method was evaluated in a moving phantom experiment, which allows the comparison with the known motion pattern and to established criteria based on anatomical images. The usefulness of the method for the calibration of optical-flow based algorithms was also demonstrated invivo under conditions similar to thermo-ablation for the abdomen of twelve volunteers. In average over all volunteers, a resulting displacement error of 1.5 mm was obtained (largest observed error equal to 4-5 mm) using a criterion based on anatomical image similarity. A better average accuracy of 1 mm was achieved using the proposed criterion (largest observed error equal to 2 mm). In both kidney and liver, the proposed criterion was shown to provide motion field accuracy in the range of the best achievable.

Index Terms-Image registration, Motion analysis, Motion compensation, Magnetic resonance imaging I. INTRODUCTION R EAL-time motion estimation on MRI is rapidly gaining importance in interventional procedures of mobile organs (such as the heart, kidney and liver). For example, miniinvasive interventional procedures such as tissue ablation with radio-frequency (RF) and non-invasive ablation techniques such as high-intensity focused ultrasound (HIFU) show a high potential in oncology and cardiology as an alternative to classical surgery [START_REF] Cline | MR-guided focused ultrasound surgery[END_REF]. Since their objective is a necrosis Sébastien Roujol, Mario Ries, Chrit Moonen and Baudouin Denis de Senneville * are with the Laboratory for Molecular and Functional Imaging: From Physiology to Therapy, FRE 3313 CNRS/University of Bordeaux 2, 33076 Bordeaux, France, (e-mail: {s.roujol,m.ries,chrit.moonen,baudouin}@imf.u-bordeaux2.fr).

Asterisk indicates the corresponding author.

Sébastien Roujol is also with the LaBRI, UMR 5800 CNRS/University of Bordeaux 1, 33405 Talence, France.

Manuscript received XX XX, 2010; revised XX XX, 2010.

of the tissue caused by an elevated thermal dose (TD) due to sustained heating, these methods are ideally paired with a non-invasive imaging modality which can provide both target tracking information and thermal dose measurements.

The most promising candidate for this role is real-time MRthermometry based on the water proton resonance frequency (PRF) [START_REF] Ishihara | A precise and fast temperature mapping using water proton chemical shift[END_REF]. This technique provides continuous temperature mapping inside the human body as well as target tracking information by exploiting the complex nature of the MRsignal: whereas the signal magnitude M provides anatomical information such as the position and the composition of the tissue, the phase ϕ is directly proportional to the local proton resonance frequency and thus the local temperature. However, for moving organs the temperature information must be mapped to a reference position in order to obtain thermaldose measurements which require the voxel-by-voxel temporal integration of the temperature. An accurate estimation of organ displacement is mandatory since spatial mismatch of few millimeters may generate a perturbation of several dozen degrees [2] [3]. In addition, when the ablation is performed with a HIFU device, the estimation of the organ displacement is the basis for the dynamic adjustment of the focal point position to track the targeted pathological tissue, which prevents an inefficient treatment and unwanted destruction of healthy tissue [START_REF] Denis De Senneville | Real time adaptive methods for treatment of mobile organs by MRI controlled high intensity focused ultrasound[END_REF].

A number of techniques have been proposed in the past to estimate on-line organ displacement, such as navigator echoes [START_REF] De Zwart | On-line correction and visualization of motion during MRI-controlled hyperthermia[END_REF] or ultrasonic echoes [START_REF] Pernot | 3D real-time motion correction in high intensity focused ultrasound therapy[END_REF]. However, the estimated motion information is generally restricted to a translational motion information. In addition, for ultrasonic echoes based approach, the estimated motion information is restricted to knowledge outside the heated zone because of the echo perturbation induced by the temperature rise. Since modern MRI acquisition methods allow the rapid acquisition of images with an excellent tissue contrast and high spatial resolution, complex deformation can be estimated on a voxel-by-voxel basis using image registration.

A variety of image registration algorithms [START_REF] Maintz | A survey of medical image registration[END_REF] has been suggested for this purpose. Since these methods are originally developed for applications based on optical imaging, they derive displacement and deformation from anatomical (i.e. magnitude) image information [START_REF] Denis De Senneville | Real time adaptive methods for treatment of mobile organs by MRI controlled high intensity focused ultrasound[END_REF]. However, the performance of these algorithms and the assessment of the quality of the obtained motion estimation is generally measured by criteria which are also based on the image magnitude similarity [START_REF] Skerl | A protocol for evaluation of similarity measures for rigid registration[END_REF]. This is for MR-imaging problematic, since the fact that estimation and control rely on the same magnitude information renders the outcome highly dependent on the noise level of the magnitude image (which can be conveniently quantified by evaluating the signal-to-noise ratio (SN R)). Unfortunately, real-time MRI of abdominal organs is frequently hampered by low SN R values of 5-20, which frequently leads to poor image registration. This problem could be alleviated if the criterion for motion estimation assessment could be based on additional information which is independent of the information used for motion estimation.

A promising way to provide this information is to exploit the phase information of the MR-image. The authors underline that it is proposed to exploit the phase component of the complex MR-signal in the image domain, not the phase associated with the fourier representation of the anatomical image [START_REF] Fleet | Computation of component image velocity from local phase information[END_REF]. While the magnitude of an MR-image reflects the underlaying anatomy, the MR-phase ϕ in gradient recalled echo images is mainly determined by the local susceptibility distribution χ and the local magnetic field strength H:

ϕ ( r) ∝ γµ 0 (1 -σ (T ( r))) (1 + χ ( r)) H ( r) (1)
whereby the gyromagnetic ratio γ and the magnetic constant µ 0 represent material constants, r the spatial position and σ the temperature (T ) dependent screening constant of the water protons. Since the susceptibility of biological tissue varies little, it is mainly the magnetic field strength H which determines the spatial variation of ϕ. Furthermore, since modern MRsystems achieve a good homogeneity of the magnetic field strength H across the field-of-view, ϕ is a spatially slowly varying function which is differentiable within homogeneous tissue. These characteristics appear advantageous for a basis of tissue displacement detection.

Consequently, the objective of this study is to introduce a new criterion to quantify the quality of the estimated motion, based on the MR-phase image similarity (P IS). This criterion is used to assess, for a set of positions observed before the interventional procedure, the optimal configuration of the registration algorithm. Since for HIFU ablations of mobile organs, motion is predominantly caused by the respiratory or the cardiac cycle and is thus periodic, this optimal configuration can be used during the interventional procedure. To demonstrate the usefulness of this criteria, an optical flow based image registration algorithm is optimized using both a magnitude image similarity criterion (M IS) and the P IS criterion for abdominal organs subjected to respiratory motion.

In general, the optimal calibration of the registration depends on several parameters (such as target motion amplitude, noise, complexity of the deformation, image resolution, ...) and thus depends on the application (targeted organ, location of the tumor, slice orientation, ...). In this paper, an image registration algorithm derived from the Horn and Schunck algorithm was optimized using our P IS criterion in a preparative calibration step.

Since HIFU ablations of mobile organs exploit the MRimage information either for the real time monitoring of the interventional process [START_REF] Roujol | Real-time MR-thermometry and dosimetry for interventional guidance on abdominal organs[END_REF] or for direct retroactive control [START_REF] Ries | Real-time 3D target tracking in MRI guided focused ultrasound ablations in moving tissues[END_REF], the latency of the registered image stream is of great importance. For typical abdominal organ motion, the delay between the actual time of displacement and when motion information is available must not exceed 100 ms [START_REF] Ries | Real-time 3D target tracking in MRI guided focused ultrasound ablations in moving tissues[END_REF]. The implemented registration algorithm was therefore designed to be compatible with the real-time constraint and ensure short image latency. Here, the real constraint is defined as the interval of time available between two successive acquisitions, which is chosen with respect to the application.

The potential of the method was evaluated in a phantom experiment. Results of the optimal image registration were compared with gold standard positions given by an external sensor. Then, the method was demonstrated in-vivo in the abdomen of twelve volunteers under free breathing, with conditions similar to a thermo-ablation. Results of the optimal image registration were compared with manually defined gold standard positions.

II. METHOD DESCRIPTION

As shown in equation ( 1) the phase of the MR-signal represents the susceptibility distribution and the magnetic field variation in the magnet. Since both change due to organ motion, we propose to use this physical information to assess the accuracy of an estimated motion field. Two main factors, related to organ displacement, may disturb the similarity between two phase images [START_REF] Poorter | The proton resonance frequency shift method compared with molecular diffusion for quantitative measurement of two dimensional time dependent temperature distribution in phantom[END_REF]:

1) The spatial mismatch due to the displacement, which leads to a spatial mismatch of the respective magnetic susceptibility distribution χ ( r). 2) Additional phase shifts generated by a modified local demagnetisation field, which is caused by a modified magnetic susceptibility distribution [START_REF] Poorter | The proton resonance frequency shift method compared with molecular diffusion for quantitative measurement of two dimensional time dependent temperature distribution in phantom[END_REF] [12] [START_REF] Schweser | A novel approach for separation of background phase in SWI phase data utilizing the harmonic function mean value property[END_REF]. To account for these phase shifts, a precise modeling of the inhomogeneous magnetic field in-vivo is required. For that purpose, the recently suggested linear phase model approach [START_REF] Roujol | Real-time MR-thermometry and dosimetry for interventional guidance on abdominal organs[END_REF] [14] assumed a simple linear magnetic field variation between the target displacement and phase variations. Each variation of the object position leads thus to a unique phase image. If phase variations with motion have been modeled, a synthetic phase map (ϕ reco ) can be constructed. This allows to quantify registration errors by evaluating the phase similarity between any new acquisition and the corresponding synthetic one.

A 3D image acquisition would be required. However, although it is well established that MR-imaging can provide motion estimates with a high spatial resolution, it is difficult in practice to acquire on-line 3D isotropic images because of the technical limitations, spatial and temporal resolution trade-offs, and low SNR associated with fast 3D acquisition sequences. Hopefully, in practice, it is possible for respiratory induced motion to choose the imaging plane direction parallel to the principal axis of the organ displacement and thus we evaluated the technique in two dimensions.

A. Implemented image registration algorithm

Optical flow algorithms, that was first applied to medical image registration in [START_REF] Meunier | A 3-D non-linear registration algorithm for brain SPECT imaging within the Talairach reference system[END_REF], allow estimating a velocity field assuming an intensity conservation during displacement, mathematically expressed by the optical flow constraint equation (OFCE):

I x u + I y v = -I t (2)
where u and v are displacement vector components, and I x , I y , I t are the spatio-temporal partial derivatives of the pixel intensity of the image.

However, a direct estimation by minimizing the deviation from the OFCE (equation ( 2)) is an under-determined problem and thus an additional constraint is required. The algorithm proposed by Horn and Schunck introduced an additional physical constraint by assuming that the motion field is smooth in the neighborhood of estimation point (right part in equation ( 3)) [START_REF] Horn | Determining optical flow[END_REF]. For that reason, this method was found to be a good candidate to estimate organ displacements during thermoablation [START_REF] Denis De Senneville | Real time adaptive methods for treatment of mobile organs by MRI controlled high intensity focused ultrasound[END_REF]. This approach seeks a motion field minimizing:

[I x u + I y v + I t ] 2 + α 2 ∇u 2 2 + ∇v 2 2 dxdy (3)
α 2 is a weighting factor designed to link the two individual metrics (intensity variation and motion regularity): while low α 2 values allow estimation of large motion amplitude, high α 2 values increase robustness against noise or possible local intensity variations not attributed to motion. The Horn and Schunck algorithm was thus applied to anatomical (i.e. magnitude) images [START_REF] Horn | Determining optical flow[END_REF]. In order to stabilize the convergence of the algorithm, a multi-resolution scheme was used [START_REF] Pratikakis | Robust multiscale deformable registration of 3D ultrasound images[END_REF] which iterates the registration algorithm from a 4-fold downsampled image (where displacements are small and the SN R is increased by the low-pass filtering inherent to the downsampling process) step-by-step to the full image resolution.

The optical-flow algorithm, applied on magnitude images, provides a motion field with a sub-pixel precision and an interpolation was required to obtain registered phase images. Due to the 2π periodicity of the phase, the spatial transformation could not be directly applied on phase images. Although this problem could be circumvented by employing a 2D phase-unwrapping step to the phase images, this remains a computationally intensive processing step which is often unstable in areas with signal discontinuities and strong susceptibility changes frequently encountered in abdominal imaging. Therefore, we applied the estimated motion to the complex MR images, to obtain registered phase images, avoiding spatial phase wraps problems.

All computationally intensive calculations were offloaded to a dedicated graphics processing unit where a detailed description of the implementation can be found in [START_REF] Denis De Senneville | An optimised multi-baseline approach for on-line MR-temperature monitoring on commodity graphics hardware[END_REF].

B. Automatic calibration of the registration

The only required intervention of the user is to define before hyperthermia, on the reference image, a region of interest encompassing the area where the registration must be optimized (we note m the associated binary mask). A reference data set was created to sample the susceptibility perturbations with motion. For most therapeutic applications, motion is caused by the respiratory or the cardiac cycle and is thus periodic. Thus, a set of K images (K = 50 was chosen in the scope of this study) covering several motion cycles of the target with a sufficient sampling density to avoid discretization errors (5-10 images per second) was acquired. A collection of registered phase images, encoding local magnetic susceptibility variations and noted ϕ r , was built and used to assess the quality of the registration. For that purpose, all phase images were registered to a common reference position on a pixel by pixel basis using the implemented registration algorithm applied on the anatomical (i.e. magnitude) image information. Registration errors were quantified by evaluating the phase similarity between any new acquisition (ϕ r ) and a synthetic one (ϕ reco ). This synthetic phase map was computed assuming a linear relation between phase variation and motion. Since the 2π periodicity of the phase would lead to a severe bias of the similarity measure, a temporal phase unwrapping on a pixel-by-pixel basis was applied between ϕ r and ϕ reco . A temporal analysis was performed to reflect the accuracy of the registration algorithm with various amplitudes and deformations likely to be encountered. For that reason, the accuracy of each estimated motion was quantified for each image of the reference dataset, and averaged to obtain the P IS criterion, defined as:

P IS = 1 K K k=1   1 N (x,y)∈m (ϕ r (x, y, k) -ϕ reco (x, y, k)) 2   (4
) Where (x, y) denotes pixel coordinates, k the image index in the time series, and N the number of pixels set to 1 in the binary mask m. ϕ reco was evaluated for each individual pixel assuming a linear phase variation along the target displacement D as follows:

ϕ reco (x, y, k) = a(x, y).D(x, y, k) + b(x, y) (5)
where D(x, y, k) is a scalar relating the displacement amplitude and orientation along the principal axis of the estimated target motion computed as follow:

D(x, y, k) = ∆X(x, y, k).V 1 + ∆Y (x, y, k).V 2 (6)
where ∆X and ∆Y denotes horizontal and vertical components of the estimated displacement, and V = (V 1 , V 2 ) is the eigen vector associated to the highest eigen value.

µ =             . . . . . . ∆X(k -1) ∆Y (k -1) ∆X(k) ∆Y (k) . . . . . .             (7)
where ∆X(k) and ∆Y (k) denote the displacement vector averaged over the mask m.

The P IS criterion was used to determine the optimal α 2 value in a preparative calibration step. For that purpose, an exhaustive enumeration of α 2 were performed (30 α 2 values were tested between 0 and 0.75) and the value minimizing the P IS criterion was selected as the optimal parameter.

C. Criterion comparison

The P IS criterion was compared with the following criteria:

a) Magnitude image similarity (MIS): The temporal average of the mean square error between each registered magnitude image (M r ) acquired during the second step and the reference one (M ref ) was computed as follows:

M IS = 1 K K k=1   1 N (x,y)∈m (M r (x, y, k) -M ref (x, y, k)) 2   (9) b
) Gold standard error (GSE): The spatio-temporal average of the Euclidean distance between the estimated motion field and a gold standard motion information was evaluated. Since in the phantom experiment, the target undergoes a translational motion, the motion was fully characterized using a navigator echo [START_REF] De Zwart | On-line correction and visualization of motion during MRI-controlled hyperthermia[END_REF]. For the in-vivo study, validation of the alignment was based on 10 landmark points, which were manually positioned and tracked over the K images in the targeted region by a staff scientist with the precision of a pixel.

D. Experimental setup

The potential of the method to calibrate optical flow based image registration algorithms was first evaluated on a phantom experiment. Subsequently, the improved performance for the on-line estimation of organ displacements is demonstrated invivo on abdominal imaging of twelve volunteers. Dynamic MR imaging was performed on a clinical Philips Achieva 1.5 T MR-system (Philips Healthcare, Best, The Netherlands). All calculations were performed on a dual processor (3.1 GHz Penryn; four cores, INTEL Santa Clara, CA, USA) workstation with 8 GB of RAM. The GPU was a NVIDIA GTX280 card with 1 GB of DRAM. The GPU implementation was realized using CUDA.

1) Phantom study: A physiological sample with relaxation times matched to the human kidney was mounted on a motorized platform to simulate an abdominal organ moving due to respiration. The applied motion pattern consisted of a periodic sinusoidal translational displacement of 20 mm of amplitude with a period of 2 s. The object displacement was monitored by an independent measure obtained by a navigator echo (0.5 mm precision), positioned parallel to the displacement, in order to get the reference displacement for the evaluation of the GSE criterion. Dynamic MR imaging was performed with a dual shot gradient recalled echo-planar sequence, which employed the following parameters (TR=30 ms, TE=15 ms, voxel size= 2×2×5 mm 3 , FOV=256×104×5 mm 3 , echo train length=25, echo spacing=1.1 ms, flip angle=20 • , bandwidth in readout direction per pixel=1777 Hz). The study was investigated for different noise levels by adding an additional Gaussian noise to the initial complex data in a separate post processing step, in order to achieve series of images with an SNR of 5 to 15.

2) In-vivo study: Dynamic MRI was performed under free breathing conditions on the kidney of 12 healthy volunteers under real-time conditions. The singleshot gradient recalled echoplanar (EPI) sequence employed the following parameters: 3000 dynamic sagittal images acquired with an imaging framerate of 10 images/s, single slice, TR=100 ms, TE=26 ms, voxel size 2.3×3.1×6 mm 3 , FOV=300×197×6 mm 3 , echo train length=63, echo spacing=0.8 ms, flip angle=35 • , bandwidth in readout direction per pixel=2085 Hz, using a four element phased array body coil. The proposed calibration method was evaluated both in the kidney and the liver individually.

3) Statistical Analysis: The significance between the optimal α 2 values obtained with the M IS, the P IS and the GSE for a variety of noise input has been evaluated for the ex-vivo experiment using an ANOVA (Analysis of Variances) in form of a F-test with significance threshold p=0.05. If the test was found significant, additional paired t-tests were applied to the data of all pairs of criteria. A significance of p=0.05 was used and corrected with the Bonferroni method.

The same statistical study was performed to assess the significance between GSE obtained with the M IS and the P IS criteria for both phantom and in-vivo experiments.

III. RESULTS

A. Phantom study

Magnitude and phase images obtained for two different phantom positions are displayed in Fig. 1. It can be observed that, contrary to magnitude images, phase images are not only shifted with the target displacement (shown by the red dashed lines), but also prone to an additional perturbation generated by a modified local susceptibility distribution (see red arrows in Fig. 1c and1d).

Fig. 2 depicts the quality of the registration evaluated for each α 2 value with the M IS (2a) and the P IS (2b) criteria. The GSE criterion is reported in dashed lines for comparison. Although the units are not the same, curve shapes can be compared. With both criteria it can be observed that the quality of the registration turned to be very poor for extreme α 2 values (close to 0 or to 1 in this example). An optimal configuration of the implemented registration algorithm was found with both criteria: the minimum of the M IS criterion was obtained for α 2 =0.175, and the minimum of the P IS criterion was obtained for α 2 =0.475. However, while the real error of the estimated displacement was 1.16 millimeters for α 2 =0.175, this value was 0.65 millimeters for α 2 =0.475.

Motion fields estimated with α 2 =0.175 and α 2 =0.475 are reported in Fig. 3c and3d, respectively. Only the motion field estimated with an optimal α 2 value evaluated using our P IS criterion matched visually the real target motion measured with the transmission line (reported in Fig. 3a).

Those results were confirmed for all tested SN R levels. α 2 values optimizing each criteria, are displayed for all tested SNR in Fig. 4a. A statistically significant difference of the optimal α 2 values obtained with the M IS, the P IS and the GSE for a variety of noise input could be observed (p<0.05). However, only the α 2 optimized with the P IS didn't shown a significant difference with the GSE optimization (p>0.05).

As observed in Fig. 2, the M IS criterion provided lower α 2 values compared to the P IS criterion for each tested SN R. α 2 values optimizing our criterion and the GSE were similar and was found to increase from 0.4 to 0.6 with the SN R.

The error of the estimated displacement obtained with an α 2 value optimized using the P IS, the M IS and the GSE criteria for several values of SN R, are reported in Fig. 4b. While the criteria based on anatomical image similarity led to a maximal error on the estimated displacement higher than one millimeter, the proposed P IS criterion reduced this value in the range of half a millimeter. The proposed criterion provides an accurate motion field for all tested SN R levels. It can be observed that the motion field accuracy was identical when the α 2 value was optimized using the P IS or the GSE criterion.

B. In-vivo study

Over the 12 human volunteers, the motion amplitude peakto-peak obtained from the landmark points was 10 mm ± 4.5 (minimum=4, maximum=18) in the kidney and 11 mm ± 4.5 (minimum=6, maximum=18) in the liver. The SN R was evaluated to 10 ± 2.5 (minimum=7, maximum=14) in the kidney and 7 ± 3 (minimum=4, maximum=14) in the liver.

Fig. 5 shows a subset of the results obtained from the in-vivo experiments. Fig. 5a and5b show the magnitude and the phase images of the liver and the kidney in their reference position. The two manually chosen masks m on the kidney and the liver (indicated in Fig. 5a by the red and blue dashed lines, respectively) were used to restrict the calibration procedure to both organs. A second set of images show magnitude and phase at a different point of the respiratory cycle (see Fig. 5c and5d), with the contours of the reference position indicated by a dashed contour line. The maximum displacement between both images sets was in this case around 12 mm. This second set was registered to the reference position with both the M IS and the P IS calibrated registration process.

The spatial distribution of the estimated displacement amplitude is reported for the M IS and the P IS criteria in Fig. 5e and 5f, respectively. As expected, larger displacements were observed in the upper part of the liver, close to the lung. A maximal error on the estimated displacement of 8 millimeters was measured on the landmark points with the M IS in the bottom and the center of the kidney and the upper part of the liver. This maximal error decreased to 4 mm when the proposed P IS criterion was used. Another convenient way to visualize the quality of the result is to subject masks, manually set on Fig. 5c, which depicts the anatomical contours, to the same transformation. This is shown in Fig. 5g and5h respectively, where the realigned masks are overlayed as a contours over the original reference image of 5a. Note, in the two zoomed regions, the registration errors arising from the low α values of the M IS-based calibration (≈ 8 mm of mismatch were observed), which are absent when the P IS was used for the calibration.

The statistical analysis over the tested volunteers of α 2 values optimizing each criteria is reported for the kidney and the liver in Fig. 6.a and 6.b, respectively. Similarly to the exvivo experiment, the M IS criterion provided lower α 2 values compared to the P IS. The error of the estimated displacement obtained with an α 2 value optimized using the P IS, the M IS and the GSE criteria are reported in Fig. 7.a (kidney) and Fig. 7.b (liver). For both kidney and liver, a statistically significant difference of the GSE could be observed for the M IS and P IS criteria in comparison to the acquisition without image registration (p<0.05). The P IS criterion performed significantly better than the M IS (p<0.05), as confirmed in the Box-and-Whisker plots in Fig. 7, which represents the GSE evaluation over all volunteers in both kidney (7.a) and (which is in the range of errors similar to GSE optimization).

(a) (b) (c) (d) (e) (f) (g) (h)

IV. DISCUSSION

Using the implemented Horn and Schunck algorithm, a poor estimation of the displacement is performed when α 2 decreased toward 0, due to instabilities of the numerical scheme. Identically, for α 2 increasing toward infinity, tested criteria indicated a poor estimation of the displacement, as the smoothness of motion constrains the velocity amplitude estimation. The GSE criterion, as a function of α 2 , reflected these properties, as reported on Fig. 2a and2b. It depicted a global minima, corresponding to an optimal registration calibration. Only the P IS matched properly this curve (see Fig. 2b).

Using the P IS criterion, the metric used for motion estimation (magnitude signal) and the metric used to evaluate the registration accuracy (phase signal) are independent. In both ex-vivo and in-vivo studies, the optimal α 2 values were systematically lower using the M IS criterion for the calibration compared to using the P IS criterion (see Fig. 4a and6a,b). This is explained by the fact that registration of the noise on magnitude images improved the M IS criterion. Although this may not be a limitation for video encoding in which SN R is generally very high, this can be problematic in the case of fast MRI sequences which display generally a lower SN R, in particular for sequences providing a high temporal resolution. Replacing the M IS criterion for the calibration of the motion estimation by the P IS criterion, reduced the average displacement error by a factor of 2 for the ex-vivo experiment and 1.5 for in-vivo experiments, respectively. The difference resulting from the use of each criterion optimization was found statistically significant in all presented experiments. This demonstrates that the proposed P IS criterion outperforms the M IS criterion for quantifying the accuracy of a motion estimation algorithm. In addition, the proposed P IS criterion was demonstrated to result in an optimal motion estimation similar to the best achievable estimate using the GSE for all tested SN R levels in the exvivo study (see Fig. 4b) and in both kidney and liver on the in-vivo study (see Fig. 7a and7b).

Several application dependent factors require an accurate calibration of the registration. With the implemented registration model, high α 2 values increased robustness against noise or possible local intensity variations not attributed to motion, but limited estimation of strong motion amplitude. In practice, the SN R was found to vary a lot between volunteers due to different coil positionings. α 2 value optimizing both the P IS and the GSE criteria decreased for low SN R values (see Fig. 4a), since less importance had to be given to grey level intensity variations in equation (3). Identically, in the invivo study, more importance was given to the displacement field regularity constraint in equation (3) for small target motion amplitudes: while α 2 values higher than 0.7 were systematically found for displacement amplitudes lower than 6 mm (2 kidneys and 1 liver were concerned), α 2 values lower than 0.25 were required for displacement amplitudes higher than 15 mm (2 kidneys and 3 livers were concerned).

It is interesting to observe in Fig. 2b that both the P IS and the GSE, as a function of α 2 , exhibits a flat zone. A range of values for α 2 allowing an accurate registration could thus be determined. This interval was similar for both the P IS and the GSE. This opens great perspectives to adjust freely α 2 in this interval, depending on the interventional application: in the present study, we used the optimal configuration of the registration algorithm for positions observed during the calibration step. Although this should be a good solution for patient under artificial breathing or for a post-processing study, α 2 may be set to the lower bound of this range for patients under free breathing, in order to allow the registration of possible larger motion amplitudes than observed during the calibration step.

The proposed P IS criterion requires an inhomogeneous magnetic field to work well, in order to obtain a spatially differentiable phase signal, which is for interventions in the abdomen for large FOVs (>10 cm) with current clinical equipment always the case. The implemented magnetic field perturbation model assumed a simple linear magnetic field variation with organ displacement. However, although this assumption holds in general for small displacements, the precision of this simple model showed several limitations in regions displaying large susceptibility variations, such as as in the vicinity of the digestive tube or in the upper part of the liver. Those effects explains the small difference between the P IS and GSE curves in Fig. 2b and4b. In addition, the linear model could be limited in the presence of through plane motion. In this case, the measured phase signal over time will not be consistent with a given tissue, and the linear relation may become inconsistent in such area. Therefore, the P IS computation may be biased in such area. However, the spatial distribution of the linear model relevance can in practice be achieved by simply mapping the fitting error (x, y) from equation [START_REF] De Zwart | On-line correction and visualization of motion during MRI-controlled hyperthermia[END_REF]. This also allows one to remove from the mask m regions where low signal levels in conjunction with complex susceptibility variations may prevent adjustment of the linear model to the phase data.

The proposed calibration method reduced the user intervention to the determination of the mask m encompassing the targeted region in the reference image. The computation time required for one image registration was 5 ms with the used material. The calibration process required the observation of thirty α 2 values, and, for each, the registration of K = 50 images and the P IS computation with equation (4). As an indication, less than ten seconds were required for the whole calibration process. The proposed method can thus be conveniently performed just before the hyperthermia procedure to optimize an optical-flow based registration algorithm. During the intervention, the optimized registration algorithm can be integrated in a computationally effective MR-data reconstruction pipeline [START_REF] Roujol | Real-time MR-thermometry and dosimetry for interventional guidance on abdominal organs[END_REF], to provide a continuously updated target position and to feed corrections into a phased-array HIFUablator [START_REF] Ries | Real-time 3D target tracking in MRI guided focused ultrasound ablations in moving tissues[END_REF]. In the future, it should also be possible to combine beam steering for an effective volume ablation with the beam adaptation required for the motion compensation [START_REF] Mougenot | Three-dimensional spatial and temporal temperature control with mr thermometry-guided focused ultrasound (MRgHIFU)[END_REF].

V. CONCLUSION Although the autocalibration based on the M IS criterion resulted in a maximal error of 4-5 mm for the in-vivo exper-iments, the P IS criterion provided better performances that reduced this maximal error below 2 mm. The resulting mean displacement error of 1.5 mm for the M IS optimization, and 1 mm for the P IS optimization in the in-vivo experiments were significantly different. No significant difference was observed between the estimated displacement with the proposed criterion and the real target motion in both ex-vivo and in-vivo experiments. Therefore, the proposed criterion based on MRphase images, which is based on a physical parameter (i.e. the magnetic field variation with displacement) was demonstrated to allow an automatic and efficient calibration of an image registration algorithm for real-time MR-guided HIFU ablations.

Although the proposed criterion was tested with a simple registration algorithm derived from the Horn and Schunck approach, it should be possible to extend to more complex image registration algorithms such as [START_REF] Nagel | An investigation of smoothness constraints for the estimation of displacement vector fields from image sequences[END_REF] [START_REF] Cornelius | Adapting optical flow to measure object motion in reflectance and X-ray image sequences[END_REF]. Results obtained also open great perspectives for an evaluation of the method with different slice positioning and orientations, and other targeted organs such as in the heart.

a

  and b are the slope and the intercept of the simple linear regression between the registered phase value ϕ reg and the target displacement D computed as follow: a(x, y) = D(x,y).ϕreg(x,y)-D(x,y).ϕreg(x,y) D(x,y) 2 .D(x,y) 2 b(x, y) = ϕ reg (x, y) -a(x, y).D(x, y)

Fig. 1 :Fig. 2 :

 12 Fig. 1: Magnitude (a,b) and phase (c,d) images obtained for two different positions of the phantom with a SN R of 10. Note that the phase is a smooth function in space and the visible phase wraps are due to the 2π periodicity of the arctan function.

Fig. 3 :

 3 Fig. 3: (a): Gold standard motion obtained from the transmission line, (b): motion field obtained with an α 2 value optimized by the M IS criterion (α 2 = 0.175), (c): motion field obtained with an α 2 value optimized by the P IS criterion (α 2 = 0.475).

Fig. 4 :

 4 Fig. 4: (a): α 2 values optimizing each criteria, for all tested SN R, (b): Gold standard errors obtained with optimal α 2 values obtained with each criteria.

Fig. 5 :Fig. 6 :Fig. 7 :

 567 Fig. 5: Typical findings of the in-vivo study: The magnitude (a) and phase (b) image acquired at the reference position are displayed with the associated processing masks m (red and blue overlays for the kidney and the liver, respectively). The corresponding shifted/deformed magnitude and phase images taken at a different point of the respiratory cycle are shown in (c) and (d). Masks m at the reference position are added as dashed line to illustrate the displacement. The spatial distribution of the estimated displacement amplitude obtained with the M IS and the P IS are displayed in (e) and (f), respectively.As described in the text, registration results for both calibration approaches are shown as registered masks of the shifted image overlayed to the reference image (solid lines in (g) and (h)). The zoomed regions highlight the spatial mismatch between the reference image and the realigned image when a M IS based calibration is used (g), which are absent when a P IS optimization is used (h).

ACKNOWLEDGMENT

The authors would like to thank the "Ligue Nationale Contre le Cancer", the "Conseil Régional d'Aquitaine", "Diagnostic Molecular Imaging" EC-FP6-project LSHB-CT-2005-512146, the "Agence Nationale de la Recherche" (project "MRgHIFU-ALKT") and Philips Healthcare.