Baudouin Denis De Senneville

Sébastien Roujol

Silke Hey

Chrit Moonen

Mario Ries

Extended Kalman Filtering for Continuous Volumetric MR-Temperature Imaging

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

Extended Kalman Filtering for continuous volumetric MR-Temperature Imaging

Baudouin Denis de Senneville, Sébastien Roujol, Silke Hey, Chrit Moonen and Mario Ries

Abstract-Real time magnetic resonance (MR) thermometry has evolved into the method of choice for the guidance of highintensity focused ultrasound (HIFU) interventions. For this role, MR-thermometry should preferably have a high temporal and spatial resolution and allow observing the temperature over the entire targeted area and its vicinity with a high accuracy. In addition, the precision of real time MR-thermometry for therapy guidance is generally limited by the available Signal to Noise ratio (SNR) and the influence of physiological noise. MR-guided HIFU would benefit of the large coverage volumetric temperature maps, including characterization of volumetric heating trajectories as well as near-and far-field heating.

In this paper, continuous volumetric MR-temperature monitoring was obtained as follows: The targeted area was continuously scanned during the heating process by a multi-slice sequence. Measured data and a priori knowledge of 3D data derived from a forecast based on a physical model were combined using an Extended Kalman Filter (EKF). The proposed reconstruction improved the temperature measurement resolution and precision while maintaining guaranteed output accuracy.

The method was evaluated experimentally ex-vivo on a phantom, and in-vivo on a porcine kidney, using HIFU heating. On the in-vivo experiment, it allowed the reconstruction from a spatio-temporally under-sampled data set (with an update rate for each voxel of 1.143 s) to a 3D dataset covering a field of view of 142.5×285×54 mm 3 with a voxel size of 3×3×6 mm 3 and a temporal resolution of 0.127 s. The method also provided noise reduction, while having a minimal impact on accuracy and latency.

Index Terms-Magnetic Resonance Imaging, Real time systems, Motion analysis.

I. INTRODUCTION

R ECENTLY, real time magnetic resonance (MR) ther- mometry [START_REF] Denis De Senneville | Magnetic resonance temperature imaging[END_REF] [2] has evolved into the method of choice for the guidance of non-invasive interventional modalities such as high-intensity focused ultrasound (HIFU) [START_REF] Hindley | MRI guidance of focused ultrasound therapy of uterine fibroids: early results[END_REF]. Noninvasive HIFU ablation of kidney and liver tumors has been demonstrated to be feasible in a clinical environment [START_REF] Dodd | Minimally invasive treatment of malignant hepatic tumors: at the threshold of a major breakthrough[END_REF] [START_REF] Stafford | Magnetic resonance temperature imaging for focused ultrasound surgery: a review[END_REF]. The principal role of MR-thermometry for the guidance of Baudouin Denis de Senneville, Mario Ries and Chrit Moonen are with the Imaging Division, UMC Utrecht, Utrecht, Netherlands, (e-mail: {b.desenneville,m.ries,c.moonen}@umcutrecht.nl).

Baudouin Denis de Senneville, Mario Ries and Chrit Moonen are also with the Laboratory for Molecular and Functional Imaging: From Physiology to Therapy, FRE 3313 CNRS/University of Bordeaux 2, 33076 Bordeaux, France.

Baudouin Denis de Senneville is also with the "Institut de Mathématiques de Bordeaux" (IMB), UMR 5251 CNRS/Université Bordeaux 1/INRIA, F-33400 Talence, France.

Sébastien Roujol is with the Cardiovascular Division, Beth Israel Deaconess, Medical Center, Harvard Medical School, Boston, United States.

Silke Hey is with Philips Healthcare, Best, Netherlands. Manuscript received ...; revised ... HIFU ablations is the continuous monitoring of the thermal intervention for an increased patient safety. MR-thermometry can be used to provide necrosis estimates and thus to determine the therapy endpoint: The cumulative temperature allows computing on-line the deposited thermal dose which, compared to the lethal dose (taken as 43 • C during 240 minutes), provides an accurate and immediate prediction of tissue necrosis [START_REF] Sapareto | Thermal dose determination in cancer therapy[END_REF]. Thermometric information can then be used for adaptive ablation strategies, which employ feedback control of the HIFU power and dynamic modifications of the HIFU trajectory [START_REF] Mougenot | Three-dimensional spatial and temporal temperature control with MR thermometry-guided focused ultrasound (MRgHIFU)[END_REF]: Since a complete destruction of the tumor is required to assure therapeutic success, efficient ablation control strategies are required exploiting both electronic beam steering and mechanical displacements of the HIFU transducer.

The feedback control algorithm can either be based on the targeted temperature, as it is required for the use of HIFU in local drug delivery applications [START_REF] Deckers | Image-guided, non-invasive, spatiotemporal control of gene expression[END_REF], or directly on necrosis estimates as is preferable for the direct thermal destruction of tumors [START_REF] Sapareto | Thermal dose determination in cancer therapy[END_REF]. However, although it is well established that MR-imaging can provide data with a high spatial resolution, it is in practice difficult to acquire 3D isotropic image volumes of a larger field of view (FOV) in real time. In general, the precision of real time MR-thermometry for therapy guidance is limited by the available Signal to Noise ratio (SNR) and the influence of physiological motion [START_REF] Conturo | Signal to noise in phase angle reconstruction: Dynamic range extension using phase reference offsets[END_REF]. As a consequence, the inherent trade-offs in MRI between spatial resolution, volume coverage, and scan time leads for fast MRI often to a sacrifice of spatial resolution and coverage in order to improve the temporal resolution. Several strategies have therefore been proposed for reducing the scan times of a Proton Resonance Frequency (PRF) sequence, such as echo-shifted gradient echo imaging [START_REF] Moonen | A fast gradient-recalled MRI technique with increased sensitivity to dynamic susceptibility effects[END_REF], SENSE [START_REF] Pruessmann | SENSE: sensitivity encoding for fast MRI[END_REF], SMASH [START_REF] Sodickson | Simultaneous acquisition of spatial harmonics (SMASH): fast imaging with radiofrequency coil arrays[END_REF], UNFOLD [START_REF] Madore | Unaliasing by Fourierencoding the overlaps using the temporal dimension (UNFOLD), applied to cardiac imaging and fMRI[END_REF] or k-t BLAST [START_REF] Tsao | k-t BLAST and k-t SENSE: dynamic MRI with high frame rate exploiting spatio temporal correlations[END_REF]. However, each method introduces new drawbacks in terms of SNR, spatial and temporal resolution.

Initially, MRI temperature monitoring of sonications has been reported in the kidney as well as in the liver in [START_REF] Daum | In vivo demonstration of noninvasive thermal surgery of the liver and kidney using an ultrasonic phased array[END_REF] and [START_REF] Kopelman | Magnetic resonance-guided focused ultrasound surgery (MRgFUS). four ablation treatments of a single canine hepatocellular adenoma[END_REF], respectively. For both studies, the MR-temperature monitoring was achieved in a single slice at an update rate of 3.4 seconds or more. However, for direct retroactive control of the heating device, MR-thermometry should have a high temporal and spatial resolution and allow observing the temperature over the entire target area and the vicinity with a high accuracy and precision. First, both the update time and the latency of the thermal information should preferably be substancially shorter [START_REF] Ries | Real-time 3D target tracking in MRI guided focused ultrasound ablations in moving tissues[END_REF]. As a consequence, the employed processing for image reconstruction must be able to complete all computations within the time interval between two suc-cessive acquisitions and achieve a short processing latency [START_REF] Roujol | Real-time MR-thermometry and dosimetry for interventional guidance on abdominal organs[END_REF]. Second, an extended spatial coverage is a necessary prerequisite to achieve:

• Retroactive control of volumetric heating strategies, based on temperature or thermal dose measurements in a predefined volumes [START_REF] Mougenot | Three-dimensional spatial and temporal temperature control with MR thermometry-guided focused ultrasound (MRgHIFU)[END_REF]. • Accurate monitoring of collateral damages such edema induced by the low temperature diffusion close to the skin [19] [20]. Similarly, the important difference of the acoustic impedance of soft tissue and bones results in an important heating of the bones by absorption of the acoustic energy [START_REF] Conturo | Signal to noise in phase angle reconstruction: Dynamic range extension using phase reference offsets[END_REF]. A recent study reported on-line volumetric MRI thermometry in pig livers with an update rate of 400 milliseconds on a volume limited to 5 slices and a relatively high precision of the temperature estimate (the median value of the standard deviation of the temperature was lower than 2 • C) [START_REF] Quesson | Real-time volumetric MRI thermometry of focused ultrasound ablation in vivo: a feasibility study in pig liver and kidney[END_REF]. To further improve the imaging performances, recent studies investigated the possibility of using Kalman filtering to exploit the physical knowledge of the heating process [22] [23]. In [START_REF] Roujol | Robust Adaptive Extended Kalman filtering for real time MR-thermometry guided HIFU interventions[END_REF], Roujol et al. used an Extended Kalman filter (EKF) to combine measured and model predicted temperature data to achieve a measurement noise reduction. A weighting factor (the Kalman gain) allowed adjusting the confidence between the predictor model and the measured data, and was calculated iteratively based on estimates of the measurement noise and the model covariance [START_REF] Kalman | A new approach to linear filtering and prediction problems[END_REF]. The EKF employed the bio-heat transfer equation (BHTE) as the model predictor and dynamically adapts the model covariance in order to achieve measurement accuracy even in regions where the parameters of the BHTE are not exactly known. Alternatively, sparse sampling approaches in k-space have been proposed, which also exploit the physical knowledge of the heating process [START_REF] Todd | Temporally constrained reconstruction applied to MRI temperature data[END_REF] to achieve off-line an artefact free reconstruction of 3D MR-temperature maps from under-sampled k-space data using sparse sampling approaches [START_REF] Todd | Reconstruction of fully 3D high spatial and temporal resolution MR temperature maps for retrospective applications[END_REF] [START_REF] Todd | Towards real-time availability of 3-D temperature maps created with temporally constrained reconstruction[END_REF].

In the current paper, a continuous volumetric MRtemperature monitoring is obtained from spatio-temporally undersampled thermometric data as follows: The targeted area is continuously scanned during the heating process by a multislice sequence. Measured data and a priori knowledge of 3D data derived from a forecast based on a physical model are combined using an EKF. The EKF improves both the temperature measurement temporal resolution and precision while maintaining guaranteed output accuracy. The method was evaluated experimentally on both an ex-vivo phantom and an in-vivo porcine kidney heated using HIFU, and demonstrated to be compatible with a real-time implementation.

II. METHOD DESCRIPTION

The proposed approach combines temperature measurement obtained from a multi-slice acquisition with a model of temperature transfer as follows:

A. Temperature modeling using the Bio Heat Transfer Equation (BHTE) model

The bio-heat transfer equation (BHTE) was employed as the model for 3D temperature prediction, which includes the applied acoustic power P , a priori knowledge of the absorption rate α, the heat diffusion coefficient D and the perfusion value w [START_REF] Pennes | Analysis of tissue and arterial blood temperatures in the resting human forearm[END_REF].

∂ ∂t T ( r, t) = α.P ( r,t) + D.∇ 2 T ( r,t) -w.T ( r,t) (1)
where ∇ 2 is the Laplacian operator, and r = (x, y, z) is the voxel coordinate. In the implemented model, coefficients P , α, D and w were assumed to be spatially and temporally invariant. The spatial distribution of the acoustic pressure field was determined using an acoustic field simulation based on Rayleigh integration over the active transducer surface. The remaining parameters of the bio-heat transfer equation were calibrated similar to the approach proposed by Mougenot et al. [START_REF] Mougenot | MRgHIFU: Feedback temperature control with automatic deduction of BHT tissue parameters[END_REF]. The temperature modeling with the BHTE was performed using a voxel size matching the reconstructed MR measurements. For this purpose, an integration of the analytic distribution of the acoustic pressure field was individually performed over each voxel in the field of view to obtain temperature prediction matching the sampling of the MR measures.

Computation times were shortened by solving the BHTE in the Fourier domain [START_REF] Quesson | Automatic control of hyperthermic therapy based on real-time Fourier analysis of MR temperature maps[END_REF], since the problem can be transformed in Fourier space into a linear differential equation, as follows:

∂ ∂t T ( k,t) + D.k 2 + w . T ( k,t) = α. P ( k,t) (2)
where T , P denotes the Fourier transform of T and P respectively and k denotes the frequency coordinates in the Fourier domain.

B. Volumetric thermal maps calculation

The discrete solution of the BHTE represents a non-linear model for data prediction in Eq. ( 2) which was addressed using the EKF formalism [START_REF] Kalman | A new approach to linear filtering and prediction problems[END_REF]. The targeted area was continuously scanned during the heating process by a multi-slice sequence. After a new 2D MRI measurement provides a temperature update T m t , the Kalman filter computed, in a first pass, the volumetric temperature prediction V - t at time t based on the last reconstructed volumetric temperature V t-1 by applying:

V - t = f (V t-1 , u t-1 ) (3)
where f represents the BHTE model, u t-1 is the control input parameter which in our case corresponds to the HIFU delivered power at time t -1.

Subsequently, the part of the thermometric data in V - t which corresponds to the currently acquired slice position, is updated with the measured temperature map T m t in order to obtain a mixed thermometric volume V m t . Then, the a priori estimate error covariance P - t was computed as follows:

P - t = A t P t-1 A T t + Q (4)
where A t is the Jacobian matrix of partial derivatives of f at

(V t-1 , u t-1
), A T t denotes the transpose t of A t , P t-1 is the a posteriori estimate error covariance at instant t -1 and Q is the process noise covariance related to the model innacuracy.

In a second pass, the EKF combined both the measured thermometric data, which is contained in V m t and the predicted data in V - t , in order to obtain the Kalman filtered temperature V t . For this purpose, a weighting factor, the Kalman gain K t , allowed to adjust in the acquired slice, the confidence between the employed model and the measured data. The Kalman gain was calculated iteratively and is predominantly influenced by two parameters: the process noise covariance Q and the measurement noise (we note R for the measurement noise covariance that relates to the MR-thermometry precision of the measure). K t is first updated based on the new a priori estimate error covariance:

K t = P t H T HP - t H T + R -1 (5)
where H relates the model state to the measurement (here H is the identity since the measurements are directly obtained in state space).

Then, the Kalman filtered temperature V t could be computed from the weighted innovation S t (S t = V m t -HV - t ) as follows:

V t = V - t + K t S t (6)
whereby S t is only nonzero at the location of the new thermometric measurements. Finally, the a posteriori estimate error covariance P t was updated as follows:

P t = (I -K t H) P - t (7)

C. Estimation of input parameters for the Kalman filter

The measurement noise R could be determined before the heating process by evaluating the temperature variance in the targeted area on a thermometric dataset in absence of any heating (in this study we assumed R to be spatially and temporally invariant). However, the process noise covariance Q, which characterizes the accuracy of the BHTE model to predict the true temperature evolution, is a priori not known and has shown to vary over time (especially between heat up and cool down period). As shown in [START_REF] Roujol | Robust Adaptive Extended Kalman filtering for real time MR-thermometry guided HIFU interventions[END_REF], this can be achieved using an adaptive EKF (AEKF), where Q is automatically adjusted over time based on a dynamic evaluation of the model accuracy for each new measurement. Based on the assumption that the temperature noise is a white noise centered around the true temperature, it was proposed to evaluate the model accuracy at time t by computing the spatio-temporal sum of the difference between predicted and measured data. In the current study, the model accuracy t was evaluated on a voxel-byvoxel basis by computing the temperature bias over a temporal window, as follows:

t (x, y, z) = 1 N t i=t-N +1 V - i (x, y, z) -V m i (x, y, z) (8)
The following heuristic was used to find a map Q minimizing the metric t on a voxel-by-voxel basis: A maximum acceptable penalty on the temperature bias (noted threshold ) was chosen. Two states were defined for the Kalman filter:

1) A model covariance Q(x, y, z) = Q min was employed
for voxels depicting a model accuracy t (x, y, z) lower than threshold . 2) Q(x, y, z) = Q max was chosen for t (x, y, z) > threshold .

A value of 0.1 was defined for Q min in the scope of this study. However, the determination of Q max is of great importance and was defined as the maximal temperature error within the update time between two successive slice updates (which is equal to the repetition time multiplied by the applied power and by the absorption rate).

D. Experimental setup

Continuous volumetric MR-temperature imaging was performed on a clinical Philips Achieva 1.5 T MRI (Philips Healthcare, The Netherlands). Heating was achieved with a Philips Sonalleve HIFU system (Philips Healthcare, Finland), which is integrated into the patient bed of the MR-system.

1) Phantom heating study: MRI guided HIFU heating was performed on a phantom using a PRF sequence which acquires one slice placed in the coronal direction, sweeping continuously through ten positions within the desired observation area (Temporal resolution=0.17 s per position). 1000 slices were acquired using a single-shot EPI sequence (T R=17 ms, T E=11 ms, flip angle=10 • , FOV=76 × 100 × 30 mm 3 , voxel size=1.5 × 1.5 × 3 mm 3 , which was reconstructed to a voxel size of 0.5 × 0.5 × 3 mm 3 to mitigate partial volume effects). Heating was performed using 50 Watts of electrical power during a period of 30 seconds. The accuracy of the MR thermometry in the heated region was evaluated using the readings of a T-type thermocouple probe as an independent reference.

2) In-vivo heating study: MRI guided HIFU heating was performed in-vivo in the kidney of a pig under general anesthesia using a PRF sequence which acquires one slice placed in the coronal direction, sweeping continuously through nine positions within the desired observation area (Temporal resolution=1.143 s per position). Each new acquired slice was used to reconstruct dynamically the volumetric temperature over the observation area, as described in II-B. In addition, the sequence continuously acquires one spatially invariant slice placed in sagittal direction through the heated region (the sagittal slice and coronal slices were acquired in an interleaved fashion). This slice provided an additional temperature measurement with a nine fold increased temporal resolution (i.e. equal to 0.127 s) and was the basis for the calculation of a gold standard temperature with a reduced noise as follows: A temporal temperature averaging was performed on a voxelby-voxel basis over each new succession of nine maps, which increased the SNR by a factor 3. The control slice was positioned in a sagittal direction several millimeters apart from the focal point, in order to limit the signal perturbation of the sweeped slice. 6000 slices were acquired using a singleshot EPI sequence (T R=63.5 ms per slice, T E=25 ms, flip angle=35 • , FOV=142.5 × 285 × 54 mm 3 , voxel size=3 × 3 × 6 mm 3 , reconstructed voxel size=3 × 3 × 6 mm 3 ).

Heating was performed using 250 Watts of electrical power during a period of 50 seconds. Note that since the kidney was static during the experiment, no correction of motion related temperature errors was required.

3) Implementation: All calculations were performed using a dual processor, INTEL 3.1 GHz Penryn, two cores, 8 GB of RAM. To reduce computation time, the volumetric temperature reconstruction was limited to a region of 32 × 32 × 9 voxels centered on the targeted area. In the scope of this study, the effect of perfusion was neglected (w = 0). The heat diffusion evaluated to 0.1 mm 2 s -1 for both presented experiments. The absorption rate was evaluated to 0.15 and 0.02 KJ -1 for the phantom and the in-vivo study, respectively. N = 10 was defined so that the temporal window size for the evaluation of the temperature bias t covered a duration of a few seconds.

threshold was set to 1 • C for the heuristic adaptation of the model covariance.

III. RESULTS

A. Phantom heating study

Fig. 1 shows an example of temperature maps extracted from a volume reconstructed at the end of the heating process. A significant reduction in measurement noise can be visually observed between the original measured MR-temperature data (the spatio-temporally under-sampled measured MRtemperature maps were reconstructed using a sliding window) and the volumetric AEKF reconstructed temperature data. Fig. 2 shows that, using the reference measurement obtained with the thermocouple probe (indicated by the white arrow in Fig. 1d), the precision of the measurement (i.e the standard deviation of the difference between the red and the black curve) was found to be improved from an initial value of ±2.4 • C with the original data (2a), to ±1.4 • C using the proposed AEKF reconstruction (2b). It can be observed however that the accuracy occasionaly decreased to 5 • C off the true value with the AEKF reconstruction.

B. In-vivo heating study

Fig. 3 shows an example of temperature maps extracted from a volume reconstructed at the end of the heating process.

Two typical examples of the temporal temperature evolution obtained in the heated region are reported: A voxel located close to the focal point was selected (referred to as voxel #1), as well as a voxel placed further in the vicinity (referred to as voxel #2). The two voxel locations are indicated in Fig. 3.

To illustrate the influence of the Q value on the reconstructed temperature, we first assumed Q to be spatially and temporally invariant and we fixed a low value in order to give a high confidence in the BHTE model (a typical value of Q = Q min = 0.1 was chosen). In voxel #1, close to the focal point, accurate results with an improved precision could be visually observed, as shown in Fig. 4a: While the original measured temperature shows up to 5 • C of temporal variation between two successive acquisitions, it remained stable over the whole experiment using the EKF approach. The precision of the MR thermometry was quantified during the cool-down period (time interval 120-280 s) using the standard deviation between MR measurements and a fitted exponential decay curve, and was improved from an original value of 0.9 • C to 0.11 • C using the additional a priori knowledge of the system physics provided by the BHTE. However, it can be observed in Fig. 4b that results were found to be less accurate in voxel #2, which was located in the periphery of the focal point. Then, we applied the proposed heuristic for the dynamic update of Q on a voxel-by-voxel basis. The temporal evolution of the estimated temperature bias t is reported for the two examined voxels in Fig. 5. Concerning voxel #1, the time period with > threshold covers a duration of ∼ 2 seconds (which represents less than 1% of the total experiment time). It can also be observed that, in this voxel, the model is found to be accurate most part of the time, except when the sonication is stopped. The time period with > threshold was found to be equal to ∼ 30 seconds for voxel #2 (which represents ∼ 11% of the total experiment time). Although the model accuracy decreases strongly upon start of the sonication, the model was found to be accurate in the cool down period. Fig. 6 shows the impact of the process noise covariance Q adaptation on the EKF reconstructed temperature data: The temperature evolution measured in voxel #2 was improved compared to results shown in Fig. 4b.

The inserts of Fig. 4 and Fig. 6 display the effect of the temporal interpolation of the applied AEKF: Although the temperature information was updated only every 9 th acquisition cycle, a continuous estimation was available for all time points, which showed, compared to the measured data, a significantly reduced measurement noise.

The estimated temperature bias t was found to be maximal during the HIFU energy deposition period. Fig. 7 shows the spatial distribution of the maximal temperature bias t reached during the experiment for each voxel of the observed field of view. Although very high t values (> 3 • C) were observable in voxels prone to low MR-signal, t remained well below 2.5 • C in the heated region.

The efficiency of both the sliding window-based and the AEKF reconstruction is evaluated in Fig. 8, using temperature informations of the control slice as a reference. In order to account for the different slice geometries, a spatial temperature average was performed over four neighboring voxels (for the control slice as well as for the reconstructed temperature volume) to compare the temperature in voxels with an identical isotropic resolution (i.e. 6 × 6 × 6 mm 3 ). The temperature evolution is reported in a single voxel located in the heated region, which is present both in the control slice and in the reconstructed volumes. While a similar accuracy was obtained with the sliding window-based and the AEKF reconstruction (an average temporal temperature difference of about 0.1 • C was measured), the precision was improved (the standard deviation of temperature difference was found to decrease from an initial value of 1.2 • C to below 0.6 • C), as shown in the inserts of Fig. 8.

Less than 15 milliseconds were required for the complete calculation of one reconstructed 3D volume on our computational platform. This rendered the method compatible with sub-second MRI and the sampling frequency was only limited by the MR-acquisition time (up to 16 Hz was achievable with the used MR-imaging protocol). The latency (i.e the delay between the actual time of the measurement and the availability of the information) is equal to the sum of the remaining acquisition time after the passage of the k-space center (equal to half of the echo time with the used MRimaging protocol), the data transport duration (which was evaluated to be approximatively equal to 10 milliseconds with the used implementation) and the image processing time. A latency lower than 60 milliseconds was thus achievable for the obtention of one 3D temperature volume with the proposed AEKF reconstruction. 3), (b): further in the vicinity of the focal point (voxel #2 in Fig. 3). The black curve shows the measured MR-temperature data, while the red curve depicts the volumetric EKF reconstructed temperature data. Note the significant reduction in measurement noise obtained with the EKF reconstruction in (a), and the loss of accuracy in (b). Fig. 5: Temporal evolution of the estimated temperature bias t in the two voxels reported in Fig. 3): voxel #1 (black line) and voxel #2 (red line).

IV. DISCUSSION

The method extends our recent work reported in [START_REF] Roujol | Robust Adaptive Extended Kalman filtering for real time MR-thermometry guided HIFU interventions[END_REF] to reconstruct continuous volumetric temperature data from spatio-Fig. 6: Temporal temperature evolution obtained with a dynamic voxel-by-voxel basis Q adaptation in voxel #2 (arrow (2) in Fig. 3b). Note the improvement of the reconstruction at the beginning of the heating process, as compared to the scenario reported in Fig. 4b in which to Q was not dynamically adapted. temporally under-sampled 3D MR-temperature maps using an Extended Kalman Filter (EKF). The bio-heat transfer equation (BHTE) was employed as the model for 3D temperature prediction, which includes the applied acoustic power, a priori knowledge of the absorption rate, the heat diffusion coefficient and the perfusion value. The proposed method is a step toward a precise volumetric MR-temperature monitoring with a high spatial and temporal resolution. Since all measured temperature values are Kalman filtered using a predictor which models the underlying energy deposition and heat evacuation processes, the measurement noise is reduced. It is in particular important to note that this increase in thermometric precision may have a direct impact on the accuracy of the estimation of the tissue necrosis, since the calculation of the thermal dose is based on a temporal integration of the temperature on a voxel-by-voxel basis [START_REF] Sapareto | Thermal dose determination in cancer therapy[END_REF].

Since the BHTE models the physical processes of heat diffusion and absorption, it requires to be applied in 3D space in order to obtain unbiased results. However, the available acquisition time is generally too limited in interventional imaging applications to obtain full 3D thermometric coverage of the heated area, especially if physiological motion has to be resolved. As a consequence, previous approaches [START_REF] Roujol | Robust Adaptive Extended Kalman filtering for real time MR-thermometry guided HIFU interventions[END_REF], applied the BHTE to 2D or severely undersampled 3D datasets, which leads to a systematic under-estimation of the heat evacuation and neglects heat inflow from adjacent slices, in particular for low temporal imaging resolution. This led to a systematic accuracy bias introduced to the filtering process. In the presented approach, the targeted area was continuously scanned during the heating process by a multi-slice sequence and thus a 3D implementation of the BHTE model was possible, which improved the filtering performance in terms of both accuracy and precision.

For the in-vivo experiment, the proposed method allowed the reconstruction of thermometric 3D datasets with an interpolated temporal resolution of 0.127 s from a spatio-temporally under-sampled data (update time per voxel 1.143 s) in a temporally consistent way: Although neighboring voxels in slice direction are in the original dataset 0.127 s apart, all voxels in the 3D map show the temperature evolution at the same point in time.

The computation times for the reconstruction of a complete 3D volume were found to be below the acquisition time of one single slice (63.5 ms on the in-vivo experiment). In addition, compared to the previously proposed approach for 3D dynamic temperature reconstruction [START_REF] Todd | Reconstruction of fully 3D high spatial and temporal resolution MR temperature maps for retrospective applications[END_REF], only the knowledge of past temperature measurements is required, which renders the method compatible with a real-time implementation.

The intrinsic limit of the proposed approach arises from its dependence on a correct representation of the underlaying heating process by the BHTE model. The parameters of the BHTE model, which are in practice often not exactly known, and are subjected to spatio-temporal changes due to the heating process. Furthermore, although the acoustic field can be precisely modeled in homogeneous media, variations in the celerity and attenuation will generally modify the predicted field in biological tissue. As a consequence, the estimate based on Rayleigh integration was sufficient for accurate filtering directly in the focal point area (see Fig. 4a), but was found inaccurate in the periphery of the focus, leading to a systematic underestimation of the observed temperature (see the insert of Fig. 4b).

These model errors can lead to a systematic accuracy bias, which could be observed when Q was fixed deliberately to give a high confidence to the BHTE model, as reported in Fig. 4. In the proposed approach, this loss of accuracy was addressed by the dynamic model covariance adaptation of the EKF, which applies a stronger weighting to the measured data in the affected regions. For this purpose, the adjustment of the Q parameter allowed handling correctly the model accuracy in the proposed filter: While high Q values lead to measurement fidelity, low Q values lead to an emphasis on modeling. Although this reduces the achievable noise removal, it prevents the filtering to introduce a systematic bias and thus a loss of accuracy in the affected regions. In the scope of the presented work, the numerical complexity induced by the increased amount of data in 3D prevented the use of an iterative optimization scheme to find the optimal Q minimizing the filter accuracy t as it was proposed in [START_REF] Roujol | Robust Adaptive Extended Kalman filtering for real time MR-thermometry guided HIFU interventions[END_REF]. Instead, a simple heuristic was introduced which allowed an adaptation of Q on a voxel-by-voxel basis based on an a priori error estimate of the process noise covariance.

Another important aspect is the latency of the Q adaptation process: In the presented implementation, this latency was determined by the size of the temporal window employed to evaluate the filter accuracy t (which covered the past 10 seconds in the presented in-vivo study). Although a high temporal window size increases the precision of the evaluation, it also increases the latency of the filter adaptation. This latency may disturb the AEKF reconstruction by raising the two following scenarios:

• In the focal point, the estimated model accuracy can be biased when the HIFU generator was shifted from ON to OFF (the estimated model accuracy reached the maximal allowed threshold t , as shown on Fig. 1 and Fig. 5). • In the vicinity of the focal point, it typically introduces a delay in the EKF adaptation at the beginning of the HIFU sonication, as it can be observed in Fig. 6 and Fig. 8b. With the proposed approach, the high temporal resolution per volume associated with a low latency allowed for a typical non-invasive measurement of the local temperature distribution, which is of particular interest for highly perfused and/or mobile organs targets such as kidney and liver tumors [START_REF] Roujol | Real-time MR-thermometry and dosimetry for interventional guidance on abdominal organs[END_REF]. As a consequence, the current approach is a step toward an adequate spatial resolution over 3D volume with sufficient temporal resolution to cope with sophisticated tumor ablation procedures which involves high energy depositions, realtime feedback control, and offline optimized 3D trajectory. However, for this purpose, several limitations will need to be addressed in future studies:

First, future works must involve the integration of a 3D acoustic pressure field simulation accounting for media with heterogeneous acoustic properties in the beam path [START_REF] Khokhlova | Focus splitting associated with propagation of focused ultrasound through the rib cage[END_REF]. Although a Fourier-based approach could be used in the scope of the current study to speed up the on-line BHTE model resolution [START_REF] Quesson | Automatic control of hyperthermic therapy based on real-time Fourier analysis of MR temperature maps[END_REF] (see Eq. ( 2)), this technique is only valid under the assumption that no heterogeneous media are present. The solution of the BHTE in heterogeneous media requires more complex numerical solutions which would potentially render a real-time volumetric AEKF reconstruction challenging.

Second, the implemented BHTE must be improved by integrating several physiological parameters into the model: One limitation of the implemented predictor model arises from the negligence of the heat evacuation due to tissue perfusion. Future works on in-vivo tissue, especially in the liver, will require a careful assessment of the corresponding error bias and potentially to take the perfusion term in Eq. (1) into account. Similarly, an assessment of the local absorption and diffusion also needs to be achieved. The a priori knowledge of these parameters is in practice difficult and often relies either on empirical data or calibration experiments [START_REF] Mougenot | Three-dimensional spatial and temporal temperature control with MR thermometry-guided focused ultrasound (MRgHIFU)[END_REF]. Moreover, although apparent coefficients can be expected to be sufficiently homogeneous within large organs such as the liver, this assumption could break down in smaller heterogeneous organs such as the kidney, or close to organ boundaries, where both become spatially variant. An accurate assessment of these parameters will improve the performance of the volumetric reconstruction by enhancing the BHTE performances, which, in turn, will reduce the need of large scale Q adaptations. The latter needs to be limited due the unavoidable Q adjustment process latency (see the delay of the correction and aliasing effects in the insert of Fig. 4b which are not completely removed by the proposed approach).

Finally, in order to take into account dynamic changes in the energy deposition and heat evacuation during the ablation (acoustic absorption, heat diffusion, tissue perfusion), the concept of volumetric sonications has been combined with retroactive feedback control of the HIFU power. However, an important challenge for HIFU ablations in organs in the upper abdomen arises from the fact that the respiratory cycle causes a continuous displacement of the target area. Therefore, the focal point of the HIFU-system must be adaptively repositioned as the organ moves with respect to the external transducer in order to avoid damaging healthy tissue and to limit acoustic energy losses. It has been demonstrated that the volumetric heating procedure can be combined with HIFU motion correction by compensating the tissue displacement for each short sonication period [START_REF] Denis De Senneville | Real time adaptive methods for treatment of mobile organs by MRI controlled high intensity focused ultrasound[END_REF]. As a consequence, future work will also involve the integration of the correction of motion artifacts in thermal maps [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] [START_REF] Denis De Senneville | MRguided thermotherapy of abdominal organs using a robust PCA-based motion descriptor[END_REF], as well as the integration into the BHTE of an update of the focal point position combined with a trajectory optimized by the feedback control strategy.

V. CONCLUSION

The proposed approach combined temperature measurements obtained from a spatio-temporally under-sampled multislice acquisition with a bio-heat transfer equation for the reconstruction of 4D MR-thermometry. Compared to the approach proposed previously [START_REF] Todd | Reconstruction of fully 3D high spatial and temporal resolution MR temperature maps for retrospective applications[END_REF], which is, in its current form, intrinsically limited to post-processing, the proposed EKF reconstruction is compatible with real-time processing and has been implemented so that the computation times for the reconstruction of a complete 3D volume were found to be below the acquisition time of one single slice, associated with a short latency.

The resulting increased spatial coverage and temporal resolution may facilitate real-time feedback control strategies, which can be based on either the temperature or the thermal dose, depending on the particular application scenario. For HIFU applications, this allows maximizing the energy deposition in the targeted area while reducing undesired heating in the near-field. Hence, the described method is a step towards effective tumor ablations while inflicting minimal damage on the healthy surrounding tissue.

Fig. 1 :

 1 Fig. 1: MR-Thermometry results obtained on the phantom heating experiment with HIFU. Temperature maps are displayed along the three directions of the space, across the center of the focal point, after 30 seconds of heating. (a), (b) and (c) show the original spatio-temporally under-sampled measured MR-temperature maps. (d), (e) and (f) show the AEKF reconstructed temperature data. (a-d), (b-e) and (c-g) display coronal, sagittal and transversal directions, respectively. The read (R), phase (P) and slice (S) directions are reported on the bottom left of each image.

Fig. 2 :

 2 Fig. 2: Comparison between the PRF-based MR-thermometry (black line) and the reference temperature measured with the thermocouple probe (red line). (a): with the original spatio-temporally under-sampled measured MR-temperature, (b): with the AEKF reconstructed temperature data.

Fig. 3 :

 3 Fig. 3: MR-Thermometry results obtained on a pig kidney during HIFU ablation. The acquisition timing diagram is detailed in (a): The coronal slice, which sweeps continuously through nine positions, is indicated in blue, while the static sagittal control slice is reported in red. Temperature maps are displayed along the three directions of the space, across the center of the focal point, after 50 seconds of HIFU heating of a porcine kidney: (b), (c) and (d) show the original spatiotemporally under-sampled measured MR-temperature maps, (e), (f) and (g) show the AEKF reconstructed temperature data. The position of the control slice is shown in (e).

Fig. 4 :

 4 Fig. 4: Temporal temperature evolution obtained with a fixed Q value (we choose Q = 0.1 in order to give a high confidence in the BHTE model). (a): in a single voxel located close to the focal point area (voxel #1 in Fig. 3), (b): further in the vicinity of the focal point (voxel #2 in Fig. 3). The black curve shows the measured MR-temperature data, while the red curve depicts the volumetric EKF reconstructed temperature data. Note the significant reduction in measurement noise obtained with the EKF reconstruction in (a), and the loss of accuracy in (b).

Fig. 7 :Fig. 8 :

 78 Fig. 7: Spatial distribution of the maximal temperature bias t reached in each voxel during the porcine kidney experiment. (a), (b) and (c) display coronal, sagittal and transversal directions, respectively.

ACKNOWLEDGMENT

The author would like to thank the European Research Council (project ERC-2010-AdG-20100317, Sound Pharma) and the Foundation InNaBioSanté (project ULTRAFITT) for funding.


M Roujol
email: s.roujol@imf.u-bordeaux2.fr

Denis De Senneville

Mrs Benois-Pineau

M Quesson
email: b.quesson@imf.u-bordeaux2.fr

Robust real time constrained estimation of respiratory motion for interventional MRI on mobile organs.

Sébastien Roujol * , Jenny Benois-Pineau, Member, IEEE, Baudouin Denis de Senneville, Mario Ries, Bruno Quesson, Chrit Moonen

Abstract-Real time magnetic resonance (MR) imaging is a promising tool for image-guided interventions. For applications such as thermotherapy on moving organs, a fine image-based compensation of motion is required in real time to allow quantitative analysis, retro-control of the interventional device, or determination of the therapy endpoint. Since interventional procedures are usually restricted to a part of the organ/tissue under study, reduced FOV imaging represents a promising way to improve spatial and / or temporal resolution. However, it introduces new challenges for the target motion estimation since structures near the target may appear transiently due to the respiratory motion and the limited spatial coverage.

In this paper, a new image based motion estimation method is proposed combining a global motion estimation with a novel optical flow approach extending the initial Horn & Schunck (H&S) method by an additional regularization term. This term integrates the displacement of physiological landmarks, which are obtained in a preparation step by pattern matching into the variational formulation of the optical flow problem. A smooth regulation of the constraint point influences is achieved using a spatial weighting function. The method was compared to the same registration pipeline employing the H&S approach. A first evaluation was performed on synthetic dataset where the accuracy of the motion estimated with the H&S method was improved by a factor of 2 using the proposed approach. An in vivo study was then realized on both the heart and the kidney of twelve volunteers. Compared to the H&S approach, a significant improvement (p<0.05) of the DICE similarity criterion computed between the reference and the registered organ positions was achieved.

Index Terms-Image registration, Motion analysis, Biomedical image processing, Magnetic resonance imaging.

I. INTRODUCTION

R EAL time MR imaging proved to be a promising candi- date for guiding non-and mini-invasive surgical interventions [START_REF] Edelman | Clinical magnetic resonance imaging[END_REF]. The additional use of quantitative MR-measurements of tissue specific properties such as MR-relaxation times, or the local temperature, allows monitoring of the therapeutic progress and helps to determine the endpoint of the intervention. However, quantitative measurements require a motion correction to enable a direct comparison on a voxelby-voxel basis between the dynamically acquired images. In the particular case of MR-guided thermal ablations, real time MR-thermometry is used to compute temperature maps inside the human body [START_REF] Rieke | Mr thermometry[END_REF]. The temperature evolution can be used as an accurate and immediate prediction of tissue necrosis [START_REF] Sapareto | Thermal dose determination in cancer therapy[END_REF] and thus requires all temperature maps registered to a common position. Moreover, the whole registration process has to be realized within the interval of two MR-acquisitions to ensure the real time condition.

Typically, interventional procedures are usually restricted to a part of the organ/tissue under study making reduced field of view (FOV) imaging desirable. This would allow improvement of the spatial and / or temporal resolution in order to decrease partial volume effects (undesirable for quantitative analysis) and to increase imaging framerate (required to observe rapid phenomena). Several strategies have been proposed toward this direction by using saturation slabs [START_REF] Haase | Localization of unaffected spins in NMR imaging and spectroscopy (LOCUS spectroscopy)[END_REF], outer volume suppression [START_REF] Wilm | Reduced field-of-view MRI using outer volume suppression for spinal cord diffusion imaging[END_REF] or interactive reduced FOV imaging [START_REF] Schaeffter | Interactive reduced FOV imaging for projection reconstruction and spiral acquisition[END_REF]. However, a reduced FOV may introduce new challenges for the target motion estimation since structures may appear transiently due to the respiratory motion and the limited spatial coverage.

A variety of motion estimation algorithms [START_REF] Maintz | A survey of medical image registration[END_REF], [START_REF] Makela | A review of cardiac image registration methods[END_REF] have been suggested in the field of medical imaging. Optical flow algorithms [START_REF] Barron | Performance of optical flow techniques[END_REF] have been proposed for motion estimation on abdominal organs for MR-guided laser ablation [START_REF] Zientara | MRI monitoring of laser ablation using optical flow[END_REF]. Recently, an approach using an initial global motion estimation followed by an optical flow algorithm was developed for real time MR-thermometry in abdominal organs [START_REF] Roujol | Real-time MR-thermometry and dosimetry for interventional guidance on abdominal organs[END_REF]. Optical flow algorithms allow estimating a velocity field assuming an intensity conservation during displacement, mathematically expressed by the optical flow equation (OFE):

I x u + I y v = -I t . (1)
Here u and v are the displacement vector components, and I x,y,t are the spatio-temporal partial derivatives of the image pixel intensity. However, a direct estimation by minimizing the deviation from OFE (equation 1) is an under-determined problem and an additional constraint is required. The method proposed by Horn and Schunck (referred to as H&S in this paper) introduces additional physical constraints enforcing the smoothness of the motion field [START_REF] Horn | Determining optical flow[END_REF]. They seek u(x, y) and v(x, y) minimizing:

E(u, v) = [I x u + I y v + I t ] 2 + α 2 ∇u 2 2 + ∇v 2 2 dxdy, (2)
where α 2 is a weighting factor designed to link both intensity variation and motion field regularity and ∇u and ∇v are the spatial gradients of u(x, y) and v(x, y), respectively. However, when reduced FOV imaging is employed, the assumption of energy conservation may be violated due to the potential presence of transient structures. In such conditions, the Horn & Schunck algorithm is expected to fail locally (near the transient structure) to recover the correct motion. This is known as one of the ill-posed problems in motion estimation, namely the occlusion problem [START_REF] Horn | Determining optical flow[END_REF]. In detail, if a structure is only present in the reference image, the algorithm might try to match the signal from this structure (in the reference image) with a different structure (in the image to be registered). To improve the robustness of the algorithm against this effect, one can increase the value of α 2 to increase the smoothness constraint of the motion field. However, this will reduce at the same time the ability of the algorithm to provide reliable motion estimation in presence of complex deformation.

Therefore, the first step of this study was to hypothesize that a combination of intensity and landmark registration could then be used to improve the robustness of the registration against transient structures. Such combinations have been proposed in the past such as in [START_REF] Cachier | Multisubject non-rigid registration of brain mri using intensity and geometric features[END_REF] where a large number of landmarks were employed for brain registration. This method was further extended to deal with a smaller number of points (manually defined) for interactive registration of medical images [START_REF] Azar | An interactive hybrid non-rigid registration framework for 3d medical images[END_REF]. In these approaches, both intensity and landmark metrics were minimized sequentially inside an iterative minimization process. Several unified minimization framework were also proposed. In [START_REF] Loncaric | Point constrained optical flow for LV motion detection[END_REF] the algorithm was designed to obtain a final motion field fitting each pre estimated landmark motion. In the presented application, due to the low SNR, a non negligible uncertainty of the landmark correspondence may be expected and would thus strongly affect such registration approach. Another solution proposed by Becciu et al. [START_REF] Becciu | A multi-scale feature based optic flow method for 3D cardiac motion estimation[END_REF], attempt to use tags obtained from MRtagged sequence in their registration framework for cardiac contraction assessment. Unfortunately, this sequence is generally unsuitable for interventional MRI, since images are tagged by regular lines where the signal has been removed. Recently, a variational approach, integrating segmented region motion, was proposed for large displacement estimation [START_REF] Brox | Large displacement optical flow[END_REF]. This method uses a linearized OFE deviation together with regularization terms which include the correspondence of region displacements in the image plane. Despite the interest of such an approach in general purpose video sequences, its application to MRI sequences is not straightforward due to the inherent difficulties of segmentation of frames into spatially coherent regions. Recently, a grid-based deformation model was proposed [START_REF] Sotiras | Simultaneous geometric-iconic registration[END_REF], but was also shown to be sensitive/limited by the landmark extraction process according to the authors.

In this paper, we propose a new real-time motion estimation method for MRI sequences which can operate in the case of reduced FOV imaging. Preliminary results of this approach can be found in [START_REF] Roujol | Real time constrained motion estimation for ECG-gated cardiac MRI[END_REF]. The contribution of the paper is twofold. We first define the two stage interventional protocol for a robust choice of constraint points and then formulate the constrained optical flow estimation by introducing an additional regularization term in the H&S method. In our formulation, the introduction of a smooth weighting function allows for a local control of the influence of constraint points. Furthermore, in order to ensure the real-time requirement together with a short latency, all computationally intensive calculations were off-loaded to a dedicated graphics processing unit (GPU). The proposed algorithm is referred to as constrained motion estimation (CME) in the scope of this paper. It was compared with the same registration pipeline employing the H&S optical flow approach. Algorithm evaluations were conducted on both synthetic data and cardiac & kidney MR-images of healthy volunteers under free breathing conditions.

II. MATERIAL AND METHODS

The proposed CME algorithm is a two-step procedure (see Fig. 1). The first step consists in selecting the constraint points along the boundary of the organ in the reference image of the time series. The contour is then subsampled and refined based on the position of extracted feature points. In a second step the motion is estimated for each image as follows: a global translational motion estimation is performed and used to initialize a local estimation of the displacement of constraint points. Non-physiological constraint point displacements are automatically identified and corresponding constraint points are discarded. The displacements of the constraint points are then integrated into the constrained optical flow algorithm (using the global estimated motion as preconditioning) to obtain the final motion field. A detailed description of each algorithm step is now presented in the next sections. Fig. 1: General scheme of the algorithm. Prior to the intervention, constraint points are automatically extracted from a reference image (step 1). Then, during the procedure, the motion field is estimated for each frame (step 2).

A. Step 1: Constraint point selection

Anatomical points are localized and tracked over the time in order to guide and to constrain the motion estimation of the target. To select them, anatomical structures such as organ boundaries, which remain present during the acquisition and follow the target, are suitable. For this, a region of interest (ROI) was manually set around the target of the reference image. The contour of the ROI was first extracted and then regularly sampled in space to obtain a set of N points surrounding the target. To allow a certain degree of freedom on the ROI drawing a refinement of the location of the N points was then performed by moving them a near feature point computed on the reference image. Due to the abundance of works in stereo matching and image retrieval, a large amount of feature point detectors has been tested and reported in literature [START_REF] Mikolajczyk | Scale and affine invariant interest point detectors[END_REF]. The critical point is the stability of these methods with respect to affine transformations of image plane, lightening, scale variations and noise. In the case of MR images, the noise and deformable motion are the main factors. According to the evaluation in [START_REF] Mikolajczyk | Scale and affine invariant interest point detectors[END_REF], the Harris-Stephens detector [START_REF] Harris | A combined corner and edge detector[END_REF] appeared to provide a good compromise between robustness and computation time. The feature point detection is based on the following response function:

R(x, y) = Det(M x,y ) -k.T r(M x,y ) 2 , (3)
with

M x,y = i,j∈S w i,j I 2 x I x I y I x I y I 2 y i,j , (4)
where (x, y) denotes the spatial coordinates, Det denotes the determinant of a matrix, T r is the trace of a matrix, w is a weighting factor (Gaussian kernel over the region S centered on (x, y)) and k is a sensitivity parameter. The value of k is generally chosen in the [0.04-0.15] interval. In our case, a value of 0.04 was empirically found suitable for our application. The response is positive in a corner region. Therefore, the feature point with the higher response in a small neighborhood (e.g. 3x3 pixels) of the contour point was selected as landmark point. The small neighborhood search was designed to prevent positioning a landmark on a different structure/organ that may depict a different motion from the targeted organ. If no feature point was present, the initial contour point was then selected.

B. Step 2: motion estimation algorithm

An optical flow based algorithm is more efficient when it is initialized near the global optimum. Therefore, to initialize it, a global motion estimation is first performed with a simple translational model. The translation parameters (horizontal and vertical) are estimated using a sign-gradient-descent with fixed step inspired by the Netravali-Robbins method [START_REF] Tekalp | Digital video processing[END_REF]. The estimation is restricted to the ROI defined in section II-A.

Then, constraint point displacements are individually estimated (two translation parameters) using the global estimated displacement as initial estimate. This estimation is restricted to a small patch centered on each constraint point intersected with the initial ROI to allow a local refinement of the global displacement. We experimentally found that a patch size of 10×10 pixels was satisfactory using a given voxel size of 2-3 mm obtained with the employed MR-acquisitions.

To remove occasional non-physiological estimates, the following outliers rejection was added. The displacement vector (d x , d y ) of a constraint point was supposed to follow a bivariate Gaussian distribution with independence of d x and d y coordinates. A constraint point was automatically rejected if at least one of its displacement components violated the marginal 3-sigma rule.

The idea of the presented approach is to constrain the H&S formulation by locally estimated displacements of feature points. Hence, we propose the following extension of the Horn & Schunck formulation with an additional regularization term:

E c (u, v) = [I x u + I y v + I t ] 2 + α 2 ∇u 2 2 + ∇v 2 2 +λ 2 N i=1 ρ(d i , R) (u -u i ) 2 + (v -v i ) 2 dxdy, (5)
where (u i , v i ) are the horizontal and vertical components of the displacement estimated for the i th constraint point. λ 2 is the regularization parameter that allows balancing between the initial behavior of the H&S algorithm and the constraint influence. ρ is a distance function, defined as

ρ(d, R) = exp (-d 2 /R 2 ) (6)
where d represents the Euclidean distance between the pixel of coordinates (x, y) and the i th constraint point, R is a bandwidth parameter. To minimize E c (u, v), we used the calculus of variation and obtained the following system:

I 2 x u + I x I y v = α 2 ∇ 2 u -I x I t + λ 2 N i=1 (ρ(d i , R)u i ) I x I y u + I 2 y v = α 2 ∇ 2 v -I y I t + λ 2 N i=1 (ρ(d i , R)v i )
.

(7) Then, we used the approximation of the Laplacian, as suggested by Horn and Schunck, where ∇ 2 u = u -u, with u the mean value of u in the neighborhood (3x3 pixels) of the estimated point [START_REF] Horn | Determining optical flow[END_REF]. Therefore, the system can be rewritten as:

a 11 u + a 12 v = b 1 a 21 u + a 22 v = b 2 , ( 8
) with                                      a 11 = I 2 x + α 2 + λ 2 N i=1 ρ(d i , R) a 12 = a 21 = I x I y a 22 = I 2 y + α 2 + λ 2 N i=1 ρ(d i , R) b 1 = α 2 u -I x I t + λ 2 N i=1 ρ(d i , R)u i b 2 = α 2 v -I y I t + λ 2 N i=1 ρ(d i , R)v i . (9)
Finally, based on the Jacobi method, the system can be solved with the following iterative scheme:

u n+1 = b n 1 a 22 -a 12 b n 2 a 11 a 22 -a 12 a 21 , v n+1 = a 11 b n 2 -a 21 b n 1 a 11 a 22 -a 12 a 21 . (10)
C. Implementation A registration to a reference frame has been preferred to a concatenation of frame to frame estimations since the latter generally lead to error accumulations, especially when high framerate is applied such as in abdominal imaging (10-15 Hz). In cardiac application, since the imaging framerate is generally limited to the cardiac frequency, the respiratory motion between two successive images generally depicts a large amplitude, devaluating the benefit of a frame to frame motion estimation. For both the H&S and the CME implementation, the multi-resolution approach from [START_REF] Pratikakis | Robust multiscale deformable registration of 3D ultrasound images[END_REF] was used, which refines the motion estimation algorithm from the 3 rd level of sub-resolution to the full resolution. We off-loaded the most time consuming task, i.e. the iterative numerical scheme of the optical flow, to a dedicated GPU.

In-vivo calibration of the employed algorithms is not an easy task since it depends on the criterion to be optimized. The calibration of the H&S method (α 2 value) has to deal with contradictory effects. A small regularization of the motion field smoothness is required to enable the estimation of complex motion and to have a globally reliable estimated motion. On the other hand, a high constraint on the motion field smoothness would reduce the registration artifact induced by intrusive structures but will at the same time limit the ability of the algorithm to estimate complex motion. This may thus deteriorate the registration in the entire organ and thus the value of such estimated motion in the perspective of quantitative analysis and clinical application. Therefore, we decided to calibrate the employed algorithms in a way to maintain optimal performance in the general case of full FOV imaging (without the presence of intrusive structures). As recently shown [START_REF] Roujol | Automatic nonrigid calibration of image registration for real time mr-guided hifu ablations of mobile organs[END_REF], a reliable in-vivo calibration of the H&S algorithm was obtained for a range of α 2 values between 0.1-0.5. A plateau was generally observed for these ranges of values and its lower bound was suggested as a good way to cope with variations of the breathing pattern (such as an amplitude variation or drift of the respiration pattern). Therefore, a α 2 value of 0.1 was employed for the H&S algorithm. We empirically found a near-optimal solution for the CME calibration by employing the following parameters:

α 2 =0.1, λ 2 =0.1, N =20 and R 2 =5.
Note that a similar optimal configuration was obtained for the synthetic dataset experiment (see result part) except for the α 2 values that was higher due the lowest complexity of the synthesized motion.

The overall algorithm was implemented in C++ and evaluated on a dual processor (INTEL 3.1 GHz Penryn, two cores). The GPU implementation was based on the Compute Unified Device Architecture (CUDA) framework [START_REF]NVIDIA CUDA Compute Unified Device Architecture -Programming Guide[END_REF] using a NVIDIA GTX280 card.

D. Experimental setup

The proposed algorithm was evaluated on both synthetic and in vivo datasets:

1) Synthetic dataset experiment: Data creation: A sequence of T(= 30) images was created. To simulate respiratory motion typically encountered on mobile organs, a periodical (period=6 frames) geometric transformation composed by a 2D translation (T x ,T y ) and scaling (S x ,S y ) was synthesized (T x = {0, 0.5, 1, 1.5, 1, 0.5} pixels,T y = {0, 2.5, 5, 7.5, 5, 2.5} pixels, S x = S y = {1, 1.03, 1.06, 1.09, 1.06, 1.03}). A signal-tonoise ratio (SN R dB ) of 1.3 was chosen to simulate a realistic acquisition (typically between 0.7 and 1.3). A structure appearing transiently in the lower part of the image was added in half of the images to simulate the effect potentially encountered with reduced FOV imaging.

Quality assessment of the motion estimation: Since the real motion (D gt = (u GT , v GT )) and the estimated motion (D = (u, v)) are available for each pixel in such synthetic dataset experiment, the measures commonly reported in the optical flow community such as the endpoint error (EE) and the angular error (AE) of the flow [START_REF] Barron | Performance of optical flow techniques[END_REF], [START_REF] Baker | A database and evaluation methodology for optical flow[END_REF] were computed, with:

EE = (u -u GT ) 2 + (v -v GT ) 2 , (11)
AE = cos -1 1 + u × u GT + v × v GT √ 1 + u 2 + v 2 1 + u 2 GT + v 2 GT . (12)
In addition, to provide additional information about the smoothness of the estimated motion field, the harmonic energy of the estimated flow [START_REF] Vercauteren | Diffeomorphic demons using itks finite difference solver hierarchy[END_REF] was reported.

2) in vivo experiments: In vivo experiments were conducted on the heart and the kidney of a total of 12 healthy volunteers. The MRI scans were performed under free breathing conditions using a 1.5 Tesla scanner (Philips Achieva/Intera, Best, The Netherlands). The volunteers were positioned head first in supine position. The MR sequences employed saturation slabs to obtain structures appearing transiently in the FOV as typically encountered when zoom imaging is used (these conditions are then referred to as "zoom imaging conditions" in the scope of this paper).

In-vivo study on the heart under zoom imaging conditions: Dynamic MRI was performed on the heart of six healthy volunteers. The acquisition sequence was ECG-gated (to observe the heart in the same cardiac phase) using a five element phased array cardiac coil. Five contiguous adjacent slices were acquired per cycle (200 cycles per scan), in short axis view, at end diastolic phase. A slice tracking technique [START_REF] Ries | Three dimensional motion compensation for real-time MRI guided focused ultrasound treatment of abdomainal organs[END_REF] was used to compensate for respiratory motion in the third dimension. Blood signal reduction was obtained using saturation slabs positioned on each side of the imaging stack. The single shot EPI sequence employed the following parameters: FOV=260 × 260 mm 2 , voxel size=2.7 × 2.7 × 7 mm 3 , echo time=20 ms, repetition time=40 ms, SENSE acceleration factor=1.6 [START_REF] Pruessmann | SENSE: Sensitivity encoding for fast MRI[END_REF]. A saturation slab was positioned underneath the extreme position of the heart (corresponding to the position at maximum respiratory displacement) to simulate zoom imaging conditions.

In-vivo study on the kidney under zoom imaging conditions: 200 frames (single slice) in coronal orientation were acquired using a four element phased array body coil. A dual shot EPI sequence employed the following parameters: FOV=200 × 400 mm 2 , voxel size=2.3 × 2.3 × 6 mm 3 , echo time=26 ms, repetition time=52 ms, flip angle=35 • . Zoom imaging conditions were achieved using a saturation slab positioned on the top of the extreme position of the kidney.

Quality assessment of the motion estimation: In such conditions, the typical amplitude of both heart and kidney motions is about 8 pixels in imaging plane between two extreme images in the respiratory cycle. Since the true motion is unknown the quality assessment of the registration was analyzed by computing the DICE similarity coefficient (DSC) [START_REF] Dice | Measures of the amount of ecologic association between species[END_REF] between the position of the organ (ROI t ) in each registered frame (t) and its position in the reference frame (ROI ref ), as follows:

DSC(t) = 2(ROI ref ∩ ROI t ) ROI ref + ROI t , (13)
Each ROI was obtained by manual segmentation. A value of 1 for DSC indicates an ideal registration (perfect ROI matching).

The harmonic energy has been also computed and reported for all tested cases to assess the energy of the deformation fields.

III. RESULTS

A. Synthetic dataset experiment

Fig. 2: Registration example obtained on the synthetic dataset experiment. The reference image (a) and the image to be registered (b) are displayed together with the registered images obtained with the CME (c) and the H&S approach using a value of α 2 of 0.5 (d), 2 (e) and 6.5 (f). The time evolution of the average endpoint error (EE) inside the target (between both dashed red ellipses) is reported in (g).

Example of synthetic images is shown in Fig. 2 where the reference image (2a, containing the underneath transient structure) and an image corresponding the maximum synthesized motion (2b, 4 th image of the cycle) are displayed. The contours (interior and exterior) of the target obtained from the reference image are reported in dashed red curves. The H&S approach was employed to register the image in (2b) using different α 2 values. When using a small value (α 2 = 0.5), the registration is strongly influence by the underneath structure leading to a severe registration artifact (yellow arrows in 2d). By using an optimally calibrated value (α 2 = 2.5, see Fig. 3a), the registration artifact was reduced but still present (yellow arrows in 2e). A higher value of α 2 (= 10) is then reported in 2f where the registration artifact introduced by the transient structure appeared substantially reduced due to higher weight on the smoothness of the motion field. However, in this case a different registration artifact was observed due to the inability of the method to handle complex motion (here the scaling effect) as shown by the yellow arrows. The CME approach provided a reliable registration where the registered target perfectly matched the reference target contour. These results were confirmed in the plot of the averaged EE values over time (computed inside the target, between the two red dashed curves) in 2g. While the H&S approach periodically failed to recover the correction motion, the CME offered more robust performances over time.

In order to better characterize the potential gain and limitations of the compared methods, their performance in term of averaged error endpoint, harmonic energy and averaged angular error were precisely investigated in function of the employed parameters (Fig. 3). As previously observed in Fig. 2d, low α 2 values (3a, 3e and 3i) provided poor performance since the methods became very sensitive to the presence of the intrusive structure. On the contrary, high values limited the ability of the algorithm to estimate complex motion and also deteriorated the motion estimates as confirmed by the convergence of the harmonic energy toward a very small value. The averaged EE values obtained with an optimal α 2 calibration were O.42 for the H&S approach (α 2 = 2.5) and 0.21 for the CME method (α 2 =6.5) showing a reduction of the averaged EE by a factor of 2. The λ 2 value influence was then investigated ((3b, 3f and 3j). As expected, small λ 2 values tends to the H&S performance and high values tends to the extrapolation of the constraint point motions (and their associated uncertainty) leading in both cases to a deterioration of the registration. A good calibration of λ 2 (=0.1, blue curve) provided a significant improvement of the motion estimation (3b, 3j). The influence of the constraint point number (N ) was then evaluated (3c, 3g and 3k). Although an optimal value was reached around 20 points (3c, 3k, blue curve), the sensitivity of the CME to this parameter was limited. Finally the influence of the bandwidth (R 2 ) provided an optimal calibration for a value of 5 ((3d, (3h, (3i, blue curve).

B. in vivo experiments

Similar results were obtained in in vivo experiments in both the heart and the kidney of healthy volunteers. Registration Fig. 3: Influence of the free parameter tuning on the H&S and the CME methods. The averaged error endpoint (EE), the harmonic energy (HE) and the average angular error (AE) are reported in the first, second and third lines of the figure respectively. Note that the averaged EE and AE computation were restricted to the target area (defined between the two red dashed ellipses in Fig. 2a). The influence of the motion field weighting (α 2 ) on both the H&S and the CME approaches are displayed in the first column. Then, other CME parameters were analyzed: the constraint point term weighting (λ 2 , second column), the number of constraint points (N , third column) and the ray influence of the constraint point (R 2 , forth column).

examples are shown for both organs in Fig. 4. By comparing the reference images (4a, 4i) with images acquired at different positions in the respiratory cycle (4b, 4j), one can observe some structures appearing transiently (see yellow arrows) due to the signal cancellation obtained from the saturation bands. In the cardiac images, the area with low signal intensity in the lower part corresponds to a saturation band that may allow for the reduction of the FOV without additional fold-over artifacts. The signal of the liver (below the heart) almost disappeared, due to the displacement induced by the respiration. In the abdominal images, two perturbations were observed in the top part of the kidney: the liver, above the kidney, partially disappeared in 4j and the intensity of the upper part of the kidney depicted a high variation. In both examples, while the registration obtained with the H&S approach (4c, 4k) was severely deteriorated in the regions near the transient structures (see red arrows), a reliable registration was obtained in the totality of the organs using the CME (4d,4l). These findings were typical for the entire sequences as confirmed by the time evolutions of the DICE similarity coefficient and the harmonic energy, respectively shown in (4e,4m) and (4f,4n). Due to the respiratory cycle, the transient structures appeared periodically in the time series and the H&S approach periodically failed to recover a reliable motion estimate leading to low Dice similarity coefficient and elevated harmonic energy values. The CME clearly outperformed the H&S approach by providing a better overall registration and more stable performance.

Over the 12 volunteers, the averaged DICE similarity coefficient (central point inside the box) obtained with the H&S method has been significantly improved (p<0.05) using the proposed CME as shown in (4g and 4o. The minimal DICE similarity coefficient values were typically very low for certain frames using the H&S method, whereas the proposed CME allowed maintaining a better performance for all the frames (around 0.92 and 0.96 for the heart dataset and the kidney dataset, respectively)). As expected, higher harmonic energy values were obtained with the H&S method since the method periodically failed to recover the motion and also because a relative lower weight of the smoothness constraint term was employed since the same α 2 value (=0.1) was employed for both algorithms.

Constraint point filtering allowed the rejection of constraint points with non-physiological estimated displacement. In av- erage over the volunteers, less than 0.37 % and 2.67 % of the constraint points were rejected for the heart and the kidney dataset registration, respectively (with a maximum of two constraint points for a frame.

C. Real time benchmarking

Benchmarking was realized for each processing step for an image sequence of spatial resolution 128 × 128. The GPU based implementation allowed a significant reduction of the whole computation time. The total computation time of the proposed approach was evaluated to 22 ms (against 87.5 ms using a CPU only implementation) and was composed by: global motion estimation (5 ms), otion estimation of constraint points (10 ms), optical flow iterative scheme (7 ms) and image registration (0.01 ms). An acceleration factor of 10 was achieved for the computation time of the iterative numerical scheme of the optical flow algorithm (equation 10).

IV. DISCUSSION

A. Performance of the proposed method

On the synthetic dataset, the proposed CME outperformed the H&S approach that failed to estimate the real motion in presence of structures appearing transiently. The endpoint error of flow obtained with the optimal H&S calibration has been reduced by a factor of 2 using the proposed CME as shown in Fig. 3a.

Furthermore, the volunteer studies confirmed the in vivo feasibility of the CME in both the heart and the kidney. A reliable registration was obtained in the totality of the organ in all frames. Using the H&S method, similar performance were also achieved in frames with similar structures as in the reference image. However, a severe degradation of its performance was observed in presence of additional intrusive structures. Although this artifact could have been reduced by increasing the regularization of the motion field smoothness, this would have in the same time decrease its ability to estimate complex motion (as can be observed in Fig. 2f and 3a using a high value of α 2 ). As recently shown [START_REF] Roujol | Automatic nonrigid calibration of image registration for real time mr-guided hifu ablations of mobile organs[END_REF], such effect is generally observed with the H&S method by using α 2 values above the interval [0.1-0.5] with the employed in-vivo images.

B. Real time feasibility of the method

MR-guidance of interventional procedures relies on the instantaneous availability of the processed images. Therefore, this limits the available computation time. In addition, Denis de Senneville et al. demonstrated in [START_REF] Denis De Senneville | Real time adaptive methods for treatment of mobile organs by MRI controlled high intensity focused ultrasound[END_REF] that large latencies have to be compensated with the help of accurate motion prediction. However, the performance of the prediction algorithm increases greatly with short latencies. Recently, in the particular case of a HIFU ablation on mobile organs, it was demonstrated that a latency inferior to 100 ms was required for the adjustment of the beam position in order to ensure an energy deposition similar to a static experiment [START_REF] Ries | Latency compensation for real-time 3D HIFU beam-steering on moving targets[END_REF]. Here, the GPU implementation offered an acceleration factor of 10 for the computation of the constrained optical flow which is in agreement with the published work [START_REF] Roujol | Real-time MR-thermometry and dosimetry for interventional guidance on abdominal organs[END_REF], [START_REF] Denis De Senneville | An optimised multi-baseline approach for on-line MRtemperature monitoring on commodity graphics hardware[END_REF]. Significant higher acceleration factor would be expected using higher resolution as shown in [START_REF] Denis De Senneville | An optimised multi-baseline approach for on-line MRtemperature monitoring on commodity graphics hardware[END_REF] since it would benefit from a larger amount of data to process (higher occupations of each processor) and a lower relative overhead. Overall, the demonstrated CPU/GPU implementation allows the acceleration of the required processing time by a factor of four and thus ensures the real time conditions with a short low latency. Further reduction of the latency may be obtained by investigating the use of more complex optimization scheme aiming to improve the convergence speed of the algorithm.

C. Calibration of the CME

Contrary to the previous works using constraint points, a comprehensive formulation of the minimization problem was proposed. In the proposed approach, the confidence into the predetermined displacement of the selected constraint points can be freely adjusted with the regularization parameter λ 2 . The quality of the obtained optical flow depends on the quality of the initial constraint point vectors, the number of constraints (N ) and the bandwidth (R) of the ρ function parameters:

• For the motion estimation of the constraint points, only a translational model was considered as it was the most robust for small patch sizes surrounding constraint points. The optimal patch size in 128×128 MRI sequences was found to be 10×10 for our images. • In order to control optical flow, the constraint points have to be placed near eventual occlusion (or problematic area). The manual choice of constraint points is not realistic during an interventional procedure, and we can only encourage the staff physician to approximately trace the contour of the ROI. Hence, the subsampling has to be sufficient in order to get a good coverage of problematic area. On the other hand, a too large number of constraint points will slow down the computational process and may degrade the registration performance by leading to a quasi interpolation of the constraint displacements and their associated uncertainties. Therefore, for the demonstrated application, N =20 was found to be near-optimal. • The bandwidth R of the ρ function regulates the influence of remote points. The large bandwidth yields a quasi interpolation of constraint point displacements over the whole image. An optimal experimental value was R= √ 5.

• Outlier rejection for constraint point vectors was found particularly useful for small patch sizes where the estimation is more sensitive to out-of-plane motion, noise, etc.

D. Limitations, clinical perspectives and future works

Due to technical limitations of fast MR acquisition sequences, extensive 3D volume imaging on mobile organs is hard to achieve. The proposed technique has thus been evaluated in the 2D case. An extension of the method toward 3D motion correction may also be considered in future works and may contribute to better correct the motion in the third dimension by reducing out-of-plane motion artifacts. However, although the proposed algorithm can be easily extended to 3D, the main challenge would likely remain in the design of a reliable 3D MR-sequence. In such acquisition, the scan time increase would render the sequence more sensitive to intra-scan motion and fat/blood signal regrowth. Echo volumar imaging [START_REF] Rabrait | High temporal resolution functional mri using parallel echo volumar imaging[END_REF] may help to decrease the scan time of such 3D acquisition, however, its associated low resolution and its robustness against the latter artifacts should be carefully investigated.

A robust formulation of equation 5 using robust metrics such as Humber function or Lagrangian [START_REF] Black | A framework for the robust estimation of optical flow[END_REF], especially in the first two terms of equation 5 should be investigated in order to account for both brightness variation and smoothness violation (motion discontinuities). Also, the integration of a median filter at each iteration step of the iterative minimization scheme has to be carefully investigated since significant improvements have been achieved in this direction [START_REF] Sun | Secrets of optical flow estimation and their principles[END_REF]. In addition, the integration of additional regularization terms such as in [START_REF] Brox | Large displacement optical flow[END_REF] should be investigated. Although these approaches may improve the motion estimation quality, the balance between robustness, simplicity (with regard to the number of free parameters) and performance has to be carefully investigated in the perspective of a clinical use. The presented framework also opens great perspectives for integration of other motion information such as navigator or ultrasonic echoes.

The method shows a promising potential for clinical integration for two reasons. First, from its simplicity it only requires a small intervention of the staff physician (only for the mask drawing that only requires few seconds ∼5s) and few free parameters to be tuned a priori. Then, the employed values of these parameters were always identical for both heart and kidney studies demonstrating the non-necessity of a re-calibration of the parameters for each sequence and the robustness against the choice of the parameters.

Finally, although the feasibility of the method has been shown on healthy volunteers, its feasibility on other organs such as the liver or on patients with, for example, irregular cardiac motion remains to be investigated. In addition, the method has been tested in conditions of a non invasive procedure, its feasibility in the presence of invasive devices, such as a catheter, will have to be evaluated in future studies.

V. CONCLUSION In this paper a new regularization constraint of the energy functional of the H&S method was presented. This approach represents a flexible solution to integrate constraint point displacements into the optical flow estimation. This extension has been demonstrated to render optical flow methods well suited to accurately estimate the motion for interventional MRI on mobile organs in presence of intrusive structures. Significant improvements were achieved compared to the Horn & Schunck approach. Finally, the use of parallel processing on affordable commodity graphics hardware demonstrates the feasibility of the algorithm in real time with very short latency required for interventional procedures.

Fig. 4 :

 4 Fig. 4: Registration result obtained in vivo in the heart and the kidney of free breathing volunteers. Registration examples in the heart (volunteer #5) and in the kidney (volunteer #6) are given where the reference images (a,i), the images to be registered (b,j) and the registered image obtained with the H&S (c,k) and the CME (d,l) methods are displayed. The time evolution of both the DICE similarity coefficient (DSC) and the harmonic energy (HE) are shown for each registration example in (e,m) and (f,n), respectively. Finally, the DSC (g,o) and HE (h,p) are reported for each volunteer over time as Box and Whisker plot where the minimum (lower point), maximum (upper point), average (point inside the box) and standard deviation (box height) values are shown. While the H&S method was disturbed by the presence of transient structures, the CME provided robust and reliable registration performance for all the frames.

This work was supported in part by Ligue Nationale Contre le Cancer, Conseil Régional d'Aquitaine, Diagnostic Molecular Imaging, Agence National de Recherche, Fondation InNaBioSanté and Philips Medical System. , Quesson and Moonen are with the Laboratory for Molecular and Functional Imaging: from Physiology to Therapy FRE 3313