Baudouin Denis De Senneville

Abdallah El Hamidi

Chrit Moonen

A direct PCA-based approach for real-time description of physiological organ deformations

Keywords:

Dynamic MR-imaging can provide functional and positional information in real-time, which can be conveniently used on-line to control a cancer therapy, e.g. using High Intensity Focused Ultrasound or Radio Therapy. However, a precise realtime correction for motion is fundamental in abdominal organs to ensure an optimal treatment dose associated with a limited toxicity in nearby organs at risk. This paper proposes a real-time direct PCA-based technique which offers a robust approach for motion estimation of abdominal organs and allows correcting motion related artifacts. The PCA was used to detect spatio-temporal coherences of the periodic organ motion in a learning step. During the interventional procedure, physiological contributions were characterized quantitatively using a small set of parameters. A coarse-to-fine resolution scheme is proposed to improve the stability of the algorithm and afford a predictable constant latency of 80 ms.

The technique was evaluated on 12 free-breathing volunteers and provided an improved real-time description of motion related to both breathing and cardiac cycles. A reduced learning step of 10 s was sufficient without any need for patient-specific control parameters, rendering the method suitable for clinical use.

I. INTRODUCTION

R ECENT developments in rapid Magnetic Resonance Imaging (MRI), associated with fast data processing strategies, now allow acquiring functional and positional information in real-time during an interventional procedure. Dynamic MRI is thereby a promising candidate to assess an online retroactive control of the therapy. For example, real-time processing of MR-images, combined with a High Intensity Focused Ultrasound system (MR-HIFU) with rapid electronic displacement of the focal point, can be used to achieve a regional temperature control [START_REF] Stafford | Magnetic resonance temperature imaging for focused ultrasound surgery: a review[END_REF]. Similarly, the recent development of integrated MRI linear accelerators (MR-LinAc), which are designed for simultaneous irradiation and MR-imaging [START_REF] Fallone | First MR images obtained during megavoltage photon irradiation from a prototype integrated linac-MR system[END_REF], shows great potential for on-line radiotherapy guidance.

Although these new techniques appear well suited for cancer therapy in vital organs such as kidney and liver, physiological displacements induced by breathing and/or cardiac activities require a precise real-time motion management to ensure:

1) A correction of motion-induced image artifacts (in particular, MR-susceptibility variations generate apparent Baudouin Denis de Senneville and Chrit Moonen are with the Imaging Division, UMC Utrecht, Utrecht, Netherlands, (e-mail: {b.desenneville,c.moonen}@umcutrecht.nl).

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

Abdallah El Hamidi is with the "Laboratoire des Sciences de l'Ingénieur pour l'Environnement", Université de La Rochelle, France, (e-mail: aelhamid@univ-lr.fr).

Manuscript received ...; revised ... temperature perturbations in the case of MR-HIFU [START_REF] Vigen | Triggered, navigated, multi-baseline method for proton resonance frequency temperature mapping with respiratory motion[END_REF], as well as geometric image distortions in the case of both MR-HIFU and MR-LinAc [START_REF] Fallone | First MR images obtained during megavoltage photon irradiation from a prototype integrated linac-MR system[END_REF]). 2) An accurate calculation of the accumulated dose (based on MR-thermometry for thermal dose assessment with MR-HIFU [START_REF] Stafford | Magnetic resonance temperature imaging for focused ultrasound surgery: a review[END_REF] or dose simulations with MR-LinAc [START_REF] Fallone | First MR images obtained during megavoltage photon irradiation from a prototype integrated linac-MR system[END_REF]). 3) A reliable beam targeting of the pathological tissue [START_REF] Brix | Threedimensional liver motion tracking using real-time two-dimensional MRI[END_REF]. To this end, several techniques are actively developed to assess motion informations in real-time: While surrogates of physiological activities can be obtained by means of a respiratory pressure belt [START_REF] Moricawa | Feasibility of simple respiratory triggering in MR-guided interventional procedures for liver tumors under general anesthesia[END_REF] or a cardiac electrocardiogram (ECG) [START_REF] Kolandaivelu | Non-invasive assessment of tissue heating during cardiac radiofrequency ablation using MRI[END_REF], displacement measurements in the vicinity of the targeted region can be provided by navigator echos [START_REF] Köhler | Spectrally selective pencil-beam navigator for motion compensation of MR-guided high-intensity focused ultrasound therapy of abdominal organs[END_REF] or ultrasonic [START_REF] Pernot | 3D real-time motion correction in high intensity focused ultrasound therapy[END_REF] echos. More recently, fast MR-acquisition protocols allow acquiring dynamically 2D images with a respectable spatial and temporal resolution and a good contrast with a clear depiction of both targeted and healthy regions. These images can be conveniently used to estimate organ displacements during the therapy using image registration algorithms.

For this purpose, the optical flow formulation of Horn & Schunck, initially proposed in the context of motion estimation in video sequence [START_REF] Horn | Determining optical flow[END_REF], was recently shown to be well adapted for the real-time estimation of elastic organ deformations [START_REF] Denis De Senneville | Real time adaptive methods for treatment of mobile organs by MRI controlled high intensity focused ultrasound[END_REF]:

E(u, v) = Ω (I x u + I y v + I t ) 2 +α ∇u 2 2 + ∇v 2 2 dΩ
(1) where Ω is the image coordinates domain, u and v the displacement vector components, I x,y,t the spatio-temporal partial derivatives of the image pixel intensity, and α a weighting factors designed to link both intensity variation (left part of Eq. ( 1)) and motion field regularity (right part of Eq. ( 1)). Combined with a multi-resolution scheme and a fast GPU (Graphics Processing Unit) implementation, Eq. ( 1) is able to assess abdominal organ displacements with an update rate of 10 Hz for MR-images [START_REF] Roujol | Real-time MR-Thermometry and Dosimetry for interventional guidance on abdominal organs[END_REF].

However, optical flow based algorithms rely on the assumption of conservation of local intensity along the trajectory. This can be unfavorable for the clinical applications described above since intensity variations may arise from changes of MR-tissue properties during the intervention [START_REF] Rieke | MR thermometry[END_REF]. Moreover, several strategies for correcting motion artifacts necessitate the quantitative and real-time characterization of individual physiological contributions such as breathing and cardiac activities [3] [11]. Therefore, it has been proposed in [START_REF] Denis De Senneville | MRguided thermotherapy of abdominal organs using a robust PCA-based motion descriptor[END_REF] to analyze the spatio-temporal coherence of the estimated displacement using a Principal Component Analysis (PCA), in order to discard, during the intervention, in real-time, in-coherent motion patterns as follows: A reduced parameterized flow model, initially computed during a preparative learning step covering several breathing cycles, is characterized using a PCA, which provides an orthonormal basis depicting the underlying characteristic patterns of the motion. We denote by K the spatial PCA compact, positive and self adjoint operator, on W (Ω) = L 2 (Ω) × L 2 (Ω), associated to the optical flow (u, v) during the learning step. The eigenvalues of K, which are nonnegative, are sorted in decreasing order and repeated a number of times equal to their multiplicity (λ 1 ≥ λ 2 ≥ • • • ≥ 0). The technique seeks the motion descriptor (d k ) k=1•••M (M is the number of employed eigenvectors) which fulfills:

T t (x, y) = u(x, y) v(x, y) = M k=1 d k . φ k (x, y) ψ k (x, y) ∀(x, y) ∈ Ω d k = u v , φ k ψ k W (Ω) (2)
where

(φ k ) k=1•••M and (ψ k ) k=1•••M
are the horizontal and vertical components of the eigenvector #k, respectively, T t (x, y) denotes the spatial transformation between the new incoming image and a reference image, and (x, y) are the pixel coordinates. For this purpose, a minimization technique was proposed in order to find the coefficients (d k ) k=1•••M that produced a flow field T t minimizing the matching error J expressed as [START_REF] Denis De Senneville | MRguided thermotherapy of abdominal organs using a robust PCA-based motion descriptor[END_REF] [14]:

J(d 1 , d 2 , • • • , d M ) = I ref -T t -1 (I) 2 (3)
where I ref is the reference and I the incoming image. However, the function J has no convexity properties and the number of its local minima is dictated by the content of the images I ref and I. The latter may induce poor estimates of the coefficients (d k ) k=1•••M and hence the PCA sensitivity phenomena (described in Appendix A) will arise if we try to consider several eigenvectors in the basis computed in the learning step for a good representation of the movement during the interventional procedure. Another drawback rises from the fact that the computation time depends on the image content:

The latency of the obtained information is thus unpredictable which limits the application of the method, especially for feedback control strategies [START_REF] Ramsay | MR thermometry in the human prostate gland at 3.0 t for transurethral ultrasound therapy[END_REF].

The current paper aims at improving the real-time estimation and quantitative characterization of physiological displacements as follows: A direct PCA-based method is proposed, which extends the original minimization method of Eq. ( 3) by formulating the determination of the coefficients (d k ) k=1•••M with the optical flow metric of Eq. [START_REF] Stafford | Magnetic resonance temperature imaging for focused ultrasound surgery: a review[END_REF]. A coarseto-fine scheme is proposed to improve the stability for organ displacements of large amplitude. The proposed method was evaluated on 12 free-breathing volunteers and its efficiency was illustrated for real-time MR-thermometry applications. It is shown that the technique provides an improved description of motion related to both breathing and cardiac activities, with a steady latency of 80 ms, during a period of 2 minutes.

II. METHOD DESCRIPTION

A. Estimation of organ displacement during the learning step

Each time a new incoming image was acquired during a learning step of 10 s, an in-house developed, freely available, software provided a 2D motion estimate using the optical flow metric of Eq. ( 1). The reader is referred to the Appendix B for a discussion of the tool, the calibration of the input parameter α and the accuracy of the estimated motion. The proposed PCA-based approach is provided as an add-on for the tool.

B. Real-time characterization of physiological organ deformations during the interventional procedure

Our proposed approach is to find, for each new incoming image during the interventional procedure, a linear PCAdecomposition which results in a motion field completing a spatial regularity constraint. For this purpose, we suggest to minimize the energy F given by:

F (d 1 , d 2 , • • • , d M ) = E M k=1 d k φ k , M k=1 d k ψ k (4)
where the functional E is given by Eq. [START_REF] Stafford | Magnetic resonance temperature imaging for focused ultrasound surgery: a review[END_REF]. By applying the Euler-Lagrange equations on a pixel-by-pixel basis, one can derive the two equations for each (x, y) ∈ Ω, as follows:

M k=1 I 2 x φ k -α 2 ∆φ k -I x I y ψ k .d k = -I x I t M k=1 I x I y φ k -I 2 y ψ k -α 2 ∆ψ k .d k = -I y I t (5)
where ∆ denotes the Laplacian operator. We introduce the following notations for simplification:

       a k 1 (x, y) = I 2 x φ k (x, y) -α 2 ∆φ k (x, y) -I x I y ψ k (x, y) a k 2 (x, y) = I x I y φ k (x, y) -I 2 y ψ k (x, y) -α 2 ∆ψ k (x, y) b 1 (x, y) = -I x (x, y)I t (x, y) b 2 (x, y) = -I y (x, y)I t (x, y)
At this point we have 2 × Card(Ω) linear equations with common unknowns (d k ) k=1•••M . The latter can be found directly through the following overdetermined linear system:

A 1 A 2    d 1 . . . d M    = B 1 B 2 (6)
with (j = {1, 2}):

(A j ) =     . . . . . . a 1 j (x, y) . . . a M j (x, y) . . . . . .     (B j ) =     . . . b j (x, y) . . .    

C. Proposed coarse-to-fine strategy

Since the Taylor approximation of the Horn & Schunck formulation of Eq. (1) holds only for small displacements, we adapted the warping theory proposed in [START_REF] Brox | High accuracy optical flow estimation based on a theory for warping[END_REF] to the resolution of Eq. ( 4) as follows: A multi-resolution scheme was performed which iterated the registration algorithm from a 8-fold downsampled image step by step to the original image resolution. An iterative refinement process was also implemented in order to update I x,y,t within each resolution level: Default results obtained for one single iteration were compared to those obtained for a full convergence of the algorithm (the variation of d 2 between two successive iterations was compared to a maximal allowed tolerance of 10 -2 in order to ensure the convergence). To achieve a global motion regularisation (left part of Eq. ( 1)), we decomposed the overall motion descriptor (noted d ′ k ) as the sum of the contribution within the currently processed resolution (d k ) and the global contribution already estimated from the lowest resolution (d g k ). The same decomposition was performed for the calculation of the laplacian of the each overall motion component (noted ∆u ′ and ∆v ′ ), so we have:

   d ′ k = d k + d g k ∆u ′ = d ′ k ∆φ k = d k ∆φ k + d g k ∆φ k ∆v ′ = d ′ k ∆ψ k = d k ∆ψ k + d g k ∆ψ k (7)
The system of Eq. ( 6) was thus rewritten as follows:

A 1 A 2    d 1 . . . d M    = B 1 B 2 + C 1 C 2 (8)
with

(C j ) j=1,2 =     . . . c j (x, y) . . .     , c 1 (x, y) = M i=1 d g i ∆φ i (x, y) c 2 (x, y) = M i=1 d g i ∆ψ i (x, y)
For each multi-resolution level, the coefficients (d k ) k=1•••M were computed using a least square approximation. The final motion descriptor was equal to the vector (d ′ k ) k=1•••M obtained at the original image resolution.

D. Separation of physiological displacement from noise contributions

The eigenvectors associated with the M largest eigenvalues were conserved using the method proposed in [START_REF] Denis De Senneville | MRguided thermotherapy of abdominal organs using a robust PCA-based motion descriptor[END_REF]: The temporal evolution of the PCA-based motion descriptor was analyzed during the learning step in the spectral domain to separate physiological motions from the noise contribution. Typical periods of the respiratory and cardiac activities are in the range of 3 -6 s and 0.5 -2 s, respectively, hence a threshold of 4 Hz was employed to discard eigenvectors coding for noise contributions. Possible values for M were iteratively enumerated until the time course of d M (i.e the coefficient in the PCA-based motion descriptor associated to the eigenvector #M ) depicts frequencies above 4 Hz exceeding 20 % of the main peak below 4 Hz.

E. Experimental validation 1) MR imaging protocol:

Dynamic MR-imaging was performed on a Philips Achieva 1.5 T (Philips Healthcare, Best, The Netherlands) under real-time conditions. An imaging frame-rate of 10 Hz was maintained on the abdomen of 12 healthy human volunteers under free-breathing conditions. The method was evaluated in 2D and the effect of through plane motion was reduced by setting the imaging plane direction parallel to the principal axis of the organ displacement. The MR-protocol was composed of a learning step of 10 s dedicated to the calculation of the eigenvectors, followed by 2 minutes devoted to mimic an interventional procedure. The MR-sequence was a single-shot gradient recalled echo-planar with the following parameters: one coronal slice, repetition time (T R)=100 ms, echo time (T E)=26 ms, bandwidth in readout direction=2085 Hz, flip angle=35 • , field of view (F OV )=256×168 mm 2 , slice thickness=6 mm, matrix=128× 84, using a four element phased array body coil. MR-images were processed with a dual processor Intel 3.1 GHz Penryn (four cores) with 16 GB of RAM. Computational intensive calculations (i.e image down&upsampling, application of a spatial transformation, filling of matrices (A j , B j , C j ) j=1,2 in Eq. ( 8)) benefited from the multi-threaded architecture.

2) Implementation of the minimization method: The proposed direct method was compared to the existing minimization technique in terms of computation time as well as precision and accuracy of the PCA representation. To clarify the benefits on the final results, identical values for M were employed for both the minimization and the direct methods (see section II-D). A Marquardt-Levenberg algorithm [START_REF] Marquardt | An algorithm for least squares estimation on nonlinear parameters[END_REF] was employed for the minimization of the functional J. The iterative process was stopped once the variation of d 2 between two successive iterations reached a user-defined tolerance of 10 -2 . A single resolution scheme was employed by default (a justification for this choice is provided in section III). The minimization method was set in optimal conditions in order to prevent the algorithm to get caught into local minima: i) The reference image I ref was chosen in the middle of the respiratory cycle in order to limit the actual amount of displacement to estimate ; ii) The motion estimation process was restricted to a manually defined region of interest (ROI), which contains the full path of the targeted organ ; iii) The motion descriptor estimated for the previously acquired image was used as a preconditioning for the new motion estimation; iv) A spatial Gaussian filter (kernel 3 × 3) was applied to the new incoming image I t . Note that tasks disclosed in items #1 and #2 were also performed for the direct method, although not required, in order to clarify the benefits on the final results.

3) Assessment of the quality of the estimated motion: For each image acquired during the interventional procedure, motion field vectors were estimated off-line using the optical flow metric of Eq. ( 1) and taken as a reference for the evaluation of the PCA representation.

To assess the quality of the estimated motion, the pixel-wise endpoint error (EE) was computed as follows:

EE = (u -u ref ) 2 + (v -v ref ) 2 (9)
where (u, v) and (u ref , v ref ) are the estimated and the reference motion estimates, respectively. 4) Illustration of the benefit of the direct PCA-based motion descriptor for real-time MR-thermometry: The MRacquisition protocol described in II-E1 is also compatible with real-time proton resonance frequency (PRF) thermometry: The PRF is measured by obtaining the phase ϕ of the MRsignal obtained with a gradient recalled sequences, which is directly proportional to the local magnetic field strength B 0 , to the local proton resonance frequency, and thus to the local temperature [START_REF] Rieke | MR thermometry[END_REF]. An estimate of temperature changes (∆T t ) inside the human body can be obtained by evaluating phase shifts between dynamically acquired phase images (ϕ t ) and a reference non-heated data set (ϕ ref ) as follows:

∆T t = (ϕ ref -ϕ t ) .k k = (γ.ω.B 0 .T E) -1 ( 10
)
where γ is the gyromagnetic ratio (≈ 42.58 MHz/Tesla), ω is the temperature coefficient (≈ 0.009 ppm/K). However, a reliable PRF-based temperature measurement on moving targets is complicated since motion related phase variation between ϕ ref and ϕ t in Eq. ( 10) may cause strong thermometry artifacts which can bias, and even mask the true temperature evolution. Several strategies have been proposed in the past to correct online these motion related errors on thermal maps and a simple technique relies on the calculation of a motion descriptor as follows: The overall phase variation can be approximated by linear phase changes of the motion descriptor (d k ) k=1•••M on a pixel-by-pixel basis [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]. In the present study, we compare thermal maps corrected using the minimization and the proposed direct methods. The temperature stability was assessed by computing the temporal temperature standard deviation on a pixel-by-pixel basis over a user defined ROI encompassing the kidney and the liver. A paired t-test was performed to assess the significance between temperature stabilities obtained with the direct and the minimization approaches (assuming a significance threshold p = 0.05).

5) Robustness in the case of local grey level intensity variations:

For the particular case of motion estimation based on magnitude images during hyperthermia, tissue modifications induced by the heating lead to a variation of the local T 1 and T 2 relaxation time and thus to local grey level intensity modifications [START_REF] Rieke | MR thermometry[END_REF]. Consequently, the condition of energy conservation in Eq. ( 1) is locally violated which may lead to incorrect motion estimates. To analyze the impact on the proposed direct method, a signal decrease, undergoing a Gaussian spatial distribution, was simulated during the complete interventional procedure in a region located around an artery. The Gaussian signal decrease had a Full Width at Half Maximum of 15 × 15 mm 2 (100% of signal loss in the central pixel) to mimic the typical in-plane lesion size achieved during the HIFU procedure reported in [START_REF] Quesson | Real-time volumetric MRI thermometry of focused ultrasound ablation in vivo: a feasibility study in pig liver and kidney[END_REF]. The resulting bias was quantified by calculating the averaged EE over time and over a region of 15 × 15 pixels (i.e 30 × 30 mm 2 with the employed pixel size) centered on the heated area, between motion estimates obtained without and with the simulated signal decrease. The significance between the averaged EE obtained with the Horn & Schunck algorithm (Eq. ( 1)), the minimization and the proposed direct method was evaluated using an ANOVA (Analysis of Variances) in the form of a F-test with a 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 threshold of p = 0.05 was used and corrected with the Bonferroni method.

III. RESULTS

Fig. 1 shows an example on a healthy volunteer (referred to as Volunteer #3) of the characterization of abdominal organ deformations using the proposed direct PCA-based motion descriptor. Fig. 1a-d report the anatomical image of the kidney and the liver in the reference position with superimposed a subset of eigenvectors estimated using the data set obtained during the 10 s of learning step. The most important contribution on the estimated motion is induced by the breathing activity, which is thus mainly featured by the eigenvector #1: A global head-foot displacement with an amplitude increased from the bottom of the kidney to the top of the liver can be observed in Fig. 1a. The associated coefficient d 1 in the PCA motion descriptor depicts a periodical temporal variation of 4-5 s (1e), which results in a main peak in the spectral domain localized around 0.2 Hz (1g). While the eigenvector #2 encodes for more local deformations in vicinity of the vertebral column and the quadratus lamborum muscle (1b), the eigenvector #4 features the movement induced by an arterial pulsations (1c). The coefficient in the PCA motion descriptor d 4 associated with the latter depicts a periodical temporal variation of 1 second (1f), which results in a primary peak in the spectral domain localized around 1 Hz (1h). Eigenvectors associated with the lowest eigenvalues encode for the noise of the motion estimation, as shown in Fig. 1d. Fig. 2 shows the averaged EE obtained in Volunteer #1 for an incremental number of eigenvectors employed in the motion characterization and a variable number of scales included in the multi-resolution scheme. While the usage of a multi-resolution scheme hampered the performance of the minimization method (see the decline between Fig. 2e and2g), the opposite phenomenon was observed using the direct method (see Fig. 2b, 2d, 2f and 2h). For this reason, a single and a four resolution scheme were employed for the remainder of the manuscript for the minimization and the direct approach, respectively. Using the minimization method, the lowest averaged EE (equal to 0.6 mm) was obtained for M = 4 (2a). For M > 4, each additional eigenvectors deteriorates the quality of the PCA representation. In contrast, using the direct method, the averaged EE was found to decrease consistently toward 0 each time a new additional eigenvector is included in the PCA decomposition (2h). Identical behaviours were observed for images acquired during the learning step (black dashed line), as well as for images acquired during the first 40 s (red line), 80 s (green line) and 120 s (blue line) of the interventional procedure. An additional error of only ≈ 0.05 mm was introduced in the PCA representation during the interventional procedure, due to the fact that the deformation is expressed using eigenvectors optimized for images acquired during the learning step. Fig. 3 confirms the superior efficiency of the proposed direct approach for all tested volunteers. Using the direct method, results obtained for a single iteration scheme within each resolution level were comparable to those obtained iterating until convergence For each pixel the amplitude of the displacement vector was computed, and each map was individually normalized between 0 and 1 for an easier visualization. Coefficients of the motion descriptor associated to eigenvectors #1 and #4 are displayed in (e) and (g), respectively. Their corresponding representation in the spectral domain are reported in (f) and (h), respectively.

(5-10 iterations were typically necessary). This demonstrates that, given a sufficient number of iterations in the direct method, the above-mentioned comments remain valid. It is also interesting to report that the superiority of the direct method was observed whatever the number of eigenvectors employed in the minimization method. Fig. 4 details a real-time benchmarking of the proposed method for the calculation of the PCA motion descriptor of an image during the interventional procedure. The required computation time logically increases with the number of eigenvectors considered in the PCA representation for both the minimization (4a) and the direct (4b) techniques. However, only the proposed direct method provided reduced and regular computation times for each image and each tested patient. It must be noted that the computation time with the minimization method was further increased in average by 60% when (h) Fig. 2: Typical findings of the efficiency of the PCA representation obtained for different number of eigenvectors employed in the motion characterization (Volunteer #1). averaged EE over time and over both the kidney and the liver is displayed using the minimization (left column) and the direct (right column) method. Results are reported for 1 (first row), (second row), 3 (third row) and 4 (fourth row) scales employed in the multi-resolution scheme.

the motion descriptor computed for the previously acquired image was not used as preconditioning for the current motion estimation. It is also interesting to report that less than half a second was mandatory at the end of the learning step for the calculation of the eigenvectors using the truncated SVD Fig. 3: Comparison of the averaged EE over time and over both the kidney and the liver obtained using the minimization (black) and the proposed direct method for each tested volunteer. The direct method was tested using a single iteration scheme within each resolution level (dark grey) and iterating until convergence (light grey), as described in section II-C. The number of eigenvectors M selected for the representation of physiological displacements is reported for each patient above the grey bars. method proposed in [START_REF] Halko | An algorithm for the principal component analysis of large data sets[END_REF]. The interventional procedure could thus be performed in direct succession to the learning step: Since images were updated with a frequency of 10 Hz with the employed MR-sequence, only 5 images had to be discarded at the begining of the interventional procedure to prevent a latency of 0.5 s. Fig. 5 illustrates the benefit of the direct PCA-based motion descriptor for real-time MR-thermometry: Thermometry artifacts caused by susceptibility variation with motion are compensated on-line (see section II-E4). Typical results obtained on Volunteer #1 are reported in Fig. 5a and5b: Compared to the minimization method, the direct technique improved the correction of motion related errors on temperature maps for all Although a moderate improvement was observed using the proposed direct method in the kidney (arrow (1)), a reduction of thermometry artifacts by up to 3 • C could be obtained in liver, especially in the upper part (arrow ( 2)) which is subjected to elastic deformations. The upper bound of the temperature precision for 75% of the liver is reported for each tested volunteer in (c). Numbers above the grey bars recall the number of employed eigenvectors M .

tested volunteers, as shown in Fig. 5c: paired t-test showed that the direct approach performed significantly better than the minimization approach in both the kidney (p < 0.001) and the liver (p = 0.002). Although a moderate averaged improvement of the temperature standard deviation of 0.15 • C was achieved in the kidney, this value reached 0.4 • C in the liver which is more prone to elastic deformations. Fig. 6 analyzes the impact of a simulated local grey level intensity variation (see the inserts of Fig. 6a), occuring during the interventional procedure. The resulting averaged EE is reported for each method and each tested volunteer in Fig. 6b. The ANOVA showed a statistically significant difference in all volunteers between the averaged EE obtained with the three tested methods (p < 0.001). It can be observed that the poorest performer was constantly the optical flow metric of Eq. ( 1). Additional paired t-tests showed that the direct approach performed significantly better than the minimization approach (p = 0.021).

IV. DISCUSSION

Compared to the existing minimization technique, the proposed direct approach brings the following advantages:

First, an inherent drawback of the minimization method rises from the fact that the computation time is dictated by the image content. This hampered the possibility to use the method for real-time monitoring applications, in which all calculations must be done whithin the interval of time available between successive acquisitions (i.e. 100 ms in the presented results). The latency of the obtained information is also an essential condition for feedback control strategies such as for example target tracking or automatic temperature control [START_REF] Ramsay | MR thermometry in the human prostate gland at 3.0 t for transurethral ultrasound therapy[END_REF]. Over the tested volunteers, a maximal number of 8 eigenvectors were selected for the characterization of physiological organ deformations using the method detailed in section II-D. An upper bound of 30 ms on the computation time required for the calculation of the PCA-based motion descriptor was thus constantly ensured using the proposed direct approach. In practice, the image acquisition duration and the required data transport time must also be taken into account to compute the latency of the obtained information. For the employed test platform, a latency of 80 ms is achievable, composed by the sum of half of the image acquisition duration (≃ 50 ms), the required data transport time (≃ 10 ms) and the image processing time (≃ 30 ms). In particular, this latency is in accordance with the requirement for real-time target tracking during an MR-guided thermal ablation, for which an upper bound of 100 ms is mandatory as shown in [START_REF] Köhler | Spectrally selective pencil-beam navigator for motion compensation of MR-guided high-intensity focused ultrasound therapy of abdominal organs[END_REF]. It is interesting to note that, for each image acquired during the learning step, the calculation of various motion descriptors with a number of eigenvectors ranged from 1 to 8 could be achieved within the interval of time available before the next image update. Therefore, no additional time consumming processings were mandatory before the interventional procedure to compute the optimal number of eigenvectors M using the method described in section II-D. Secondly, the minimization approach, which relies on an iterative process, imposes in practice several restrictions to stabilize the convergence: Possible fold-over MR-artifacts as well as mixtures of static/dynamic parts of the entire fieldof-view and/or complex motion paterns may lead to local minima. As a consequence, the motion estimation process has to be restricted to a ROI, manually set at the beginning of the intervention, which contains the full path of the targeted organ. The proposed direct approach gets rid of this requirement. In addition, using the minimization method, motion descriptor estimated for the previously acquired image had to be used as a preconditioning for the new motion estimation. Any error in the estimation process of one image could consequently have a potential impact on the motion descriptors calculated during the rest of the intervention. This is not the case anymore with the proposed direct method where all processings are done individually for each image.

Finally, the minimization method seeks the (d k ) k=1•••M optimizing a Mean Square Error (M SE) criterion between the reference and the transformed incoming image: In other words, the image matching is in this case similar to a "visual" assessment of the quality of the registration. Using the direct approach, the estimation process during the interventional procedure is consistent with the model employed during the learning step, which leads to an improved PCA representation (Fig. 2). This allows improving in turn the correction of motion related errors on functionnal images (Fig. 5) and brings an improved robustness to possible grey level intensity variations not attributed to motion (Fig. 6).

The usage of a multi-resolution scheme had an opposite impact on the performance for the minimization and for the direct methods: Since basis flows at coarse scales are subsampled versions of the eigenvectors computed at the original scale, they may deviate slightly from orthogonality and thus be inadequate to solve the minimization problem of Eq. (3) (see Fig. 2g), as reported by Black et al in [START_REF] Black | Learning parameterized models of image motion[END_REF]. This drawback is reduced using the direct approach since, in this case, the spatial motion regularization constraint prevents such numerical instabilities (see Fig. 2h). The usage of the multi-resolution scheme was vital for the direct method in order to address the fact that the Taylor approximation of the Horn&Schunck formulation of Eq. ( 1) holds only for small displacements. Using a single resolution scheme, the direct method could not cope with the maximum amplitude of the observed organ displacements (see Fig. 2a), which were substantially bigger than the pixel size (10 mm ± 4.5 and 11 mm ± 4.5 were observed in the kidney and the liver, respectively).

It is interesting to note that such a real-time characterization of physiological motion contributions opens great perspectives for the correlation of the motion descriptor with external sensors [START_REF] Moricawa | Feasibility of simple respiratory triggering in MR-guided interventional procedures for liver tumors under general anesthesia[END_REF], [START_REF] Kolandaivelu | Non-invasive assessment of tissue heating during cardiac radiofrequency ablation using MRI[END_REF], [START_REF] Feinberg | Hybrid ultrasound MRI for improved cardiac imaging and real-time respiration control[END_REF]: The efficiency of the estimated organ deformation may be assessed by evaluating the coherence of the PCA descriptors with those independent motion informations. That way, the interventional procedure may immediately be stopped once incoherent patterns are identified in order to ensure the security of the patient.

It must be pointed out that several shortcomings persist for both minimization and direct methods: First, the evaluation was conducted on healthy subjects and the potential extension must still be assessed for patients whose physiological activities are likely to be less periodic and reproducible. It must be underlined that the proposed correction is not entirely constrained to positions present in the collection, as it can also interpolate intermediate positions. However, the observation of a significant deviation of the learned motion pattern may lead to a complete re-calibration when the change is not reversible. Moreover, the proposed technique is clearly limited in case of new observed positions during the interventional process, for which the learned motion pattern is insufficient. A combination of the method with a correction adapted for spontaneous motion (such as that described 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]) should thus be investigated in future work. Secondly, the effect of through plane motion is a serious limiting factor: Although the motion trajectory of the kidney and the lower part of the liver can be approximated in first order by a linear shift, the true trajectory in the upper part of the liver is a curve in 3D space. A dynamic 3D imaging would be optimal but the MR-acquisition is currently not fast enough to achieve the temporal resolution required to avoid intra-scan motion. Two strategies may be investigated: 1) Additional information in the third dimension, such as navigator echoes, may be used in combination with adaptive slice tracking as proposed in [START_REF] Köhler | Spectrally selective pencil-beam navigator for motion compensation of MR-guided high-intensity focused ultrasound therapy of abdominal organs[END_REF], [START_REF] Feinberg | Hybrid ultrasound MRI for improved cardiac imaging and real-time respiration control[END_REF]; 2) 3D trajectories may be estimated from 2D MRI using one or several volumetric scans obtained before the intervention as shown in [START_REF] Brix | Threedimensional liver motion tracking using real-time two-dimensional MRI[END_REF], [START_REF] Arnold | 3D organ motion prediction for MR-guided high intensity focused ultrasound[END_REF], [START_REF] Stemkens | Optimizing 4D-MRI data sampling for respiratory motion analysis of pancreatic tumors[END_REF].

V. CONCLUSION

This paper proposes a real-time PCA-based method which provides an efficient quantitative description of physiological organ deformations, with an improved steady latency, during a period of 2 minutes. The effectiveness of the method was demonstrated for real-time MR-thermometry application: An improved correction of motion related artifacts was obtained with an increased robustness to local grey level intensity variations not attributed to motion. For this purpose, a reduced learning step of 10 s was mandatory and no patient-specific control parameters needed to be set, which renders the method suitable for clinical use.

APPENDIX A PCA ANALYSIS AND BASIS SENSITIVITY

In this appendix, we will denote by K ℓ (resp. K i ) the spatial PCA compact associated to the optical flow obtained during the learning step of N images (u ℓ (x, y), v ℓ (x, y)) and the interventional step (u i (x, y), v i (x, y)), respectively. Similarly, we denote by ϕ ℓ k , ψ ℓ k k≥1 (resp. ϕ i k , ψ i k k≥1 ) the orthonormal basis of the operator K ℓ (resp. K i ). The error estimates of the PCA representation during the learning (ε l ) and the interventional (ε i ) steps can be computed using:

       ε ℓ = u ℓ v ℓ - N k=1 d ℓ k ϕ ℓ k ψ ℓ k W (Ω) = k≥N +1 λ ℓ k ε i = u i v i - N k=1 d i k ϕ i k ψ i k W (Ω)
= k≥N +1 λ i k with ε ℓ -→ 0 and ε i -→ 0 as N -→ +∞. However, in our context, we would like to expand the use the learning step basis ϕ ℓ k , ψ ℓ k k≥1 to the interventional procedure. The error of the PCA representation during the interventional procedure (ε) can be expressed as follows:

ε = u i v i - N k=1 d i,ℓ k ϕ ℓ k ψ ℓ k W (Ω) (11)
d i,ℓ k = u i v i , ϕ ℓ k ψ ℓ k W (Ω)
In this case we still have an error decrease with respect to the number of terms N (ε i -→ 0 as N -→ +∞), but the estimate ε ≤ k≥N +1 λ i k does not hold true and consequently the convergence speed of the approximation ( 11) is lost. On the other hand, the sensitivity of the PCA basis ϕ ℓ k , ψ ℓ k k≥1 , due to the clustering of the corresponding eigenvalues, can lead to bad estimates of the coefficients d i,ℓ k in the following expansion:

u i v i = k≥1 d i,ℓ k ϕ ℓ k ψ ℓ k ( 12
)
More precisely, the authors showed recently in [START_REF] Hamidi | On eigenelements sensitivity for compact self-adjoint operators and applications[END_REF] the following result: For every k ≥ 1, it holds:

|λ ℓ k -λ i k | ≤ K ℓ -K i (13)
and

ϕ ℓ k ψ ℓ k - ϕ i k ψ i k ≤ K ℓ -K i min |λ ℓ k -λ ℓ k |, |λ ℓ k -λ ℓ k | (14)
where λ ℓ k (resp. λ ℓ k ) is the closest eigenvalue among all eigenvalues which are smaller (resp. larger) than λ ℓ k . Note that this result is very classical when the eigenvalue λ ℓ k is simple and we can refer the interested reader to [START_REF] Rousselet | Continuité et différentiabilité d'éléments propres: Application à l'optimisation de structures[END_REF], [START_REF] Bhatia | The Hoffman-Wielandt inequality in infinite dimensions[END_REF]. The general case of multiple eigenvalues is studied in [START_REF] Hamidi | On eigenelements sensitivity for compact self-adjoint operators and applications[END_REF], [START_REF] Davis | The rotation of eigenvectors by perturbation. iii[END_REF]. The estimate of Eq. ( 14) is optimal in the sense that equality holds true in general situations. Moreover, when the eigenvalues λ ℓ n , λ ℓ n and λ ℓ n cluster, the eigenvector ϕ ℓ k , ψ ℓ k is very sensitive and its use in the approximation of (u i , v i ) requires an acute computation of the corresponding coefficient.

In the current paper, it is shown that local minima in the minimization method employed in [START_REF] Denis De Senneville | MRguided thermotherapy of abdominal organs using a robust PCA-based motion descriptor[END_REF], [START_REF] Black | Learning parameterized models of image motion[END_REF] will induce bad estimates of the coefficients d i,ℓ k k=1•••N and hence sensitivity phenomena, described by Eq. ( 14), will arise if we try to consider several eigenvectors in the basis (ϕ ℓ k , ψ ℓ k ) k=1•••N , for a good representation of the movement.

APPENDIX B ESTIMATION OF ORGAN DISPLACEMENT USING

OPTICAL-FLOW

The RealTITracker toolbox1 provided 2D motion estimates using the optical flow metric of Eq. [START_REF] Stafford | Magnetic resonance temperature imaging for focused ultrasound surgery: a review[END_REF]. A multi-resolution scheme was employed. To ensure the convergence of the algorithm, the averaged variation of the estimated motion amplitude was compared to a maximal allowed tolerance of 10 -3 pixels. The reader is referred to [START_REF] Roujol | Automatic non-rigid calibration of image registration for real time interventional MRI of mobile organs[END_REF] for a complete analysis of the impact of the α value on the outcome of the optical flow metric: While an increased α value intrinsically improves the robustness against low SNR values, it also limits the estimation of elastic deformations. A compromise must consequently be found. In [START_REF] Roujol | Automatic non-rigid calibration of image registration for real time interventional MRI of mobile organs[END_REF], the accuracy of the motion estimates was assessed ex-vivo using gold standard displacements provided by external sensors, and in-vivo using gold standard landmark points manually positioned and tracked by a staff scientist. For the employed implementation and the used MR-acquisition sequence, it was shown that any value in the range of 0.3 and 0.5 for α provided tracking performances within the gold standard precision for SN R of 10-15. A fixed value of 0.4 was consequently employed in the current paper.

Fig. 1 :

 1 Fig.1: Typical results obtained in the abdomen of a healthy volunteer (Volunteer #3). Kidney and liver are delimited by the white dashed line in (a) and indicated by arrows (1) and (2), respectively. (a)-(d) show a subsample of eigenvectors obtained with the displacement fields estimated during the 10 s of learning step. For each pixel the amplitude of the displacement vector was computed, and each map was individually normalized between 0 and 1 for an easier visualization. Coefficients of the motion descriptor associated to eigenvectors #1 and #4 are displayed in (e) and (g), respectively. Their corresponding representation in the spectral domain are reported in (f) and (h), respectively.

  s of learning step First 40 s of the intervention First 80 s of the intervention During the complete intervention (120 s)

Fig. 4 :

 4 Fig. 4: Boxplots of the computation time required to perform the estimation of the PCA-based motion descriptor for one image: Each box plot relates the distribution of the computation time measured on each individual image over the 2 minutes of acquisition and all tested volunteers. The boxplots are reported for different number of eigenvectors conserved in the PCA representation. The results with the minimization method (a) are compared to those obtained with the proposed direct method (b). The median is shown by the central mark, the first and the third quartile are reported by the edges of the box, the whiskers extend to the most extreme time points which are not considered as outliers, and outliers are individually plotted in red.

Fig. 5 :

 5 Fig. 5: Illustration of the benefit of the direct PCA-based motion descriptor for real-time MR-thermometry: Here, thermometry artifacts caused by susceptibility variation with motion are compensated on-line. The temperature stability is analyzed by computing the temporal temperature standard deviation on a pixel-by-pixel basis over the 2 minutes of the intervention step. The obtained maps are compared when the minimization (a) and the direct (b) PCA representations are employed for the description of abdominal displacements in volunteer #1.Although a moderate improvement was observed using the proposed direct method in the kidney (arrow (1)), a reduction of thermometry artifacts by up to 3 • C could be obtained in liver, especially in the upper part (arrow (2)) which is subjected to elastic deformations. The upper bound of the temperature precision for 75% of the liver is reported for each tested volunteer in (c). Numbers above the grey bars recall the number of employed eigenvectors M .

Fig. 6 :

 6 Fig. 6: Evaluation of the impact of a local grey level intensity variation, occuring during the interventional procedure, in the estimated motion. (a): Anatomical image obtained for volunteer #3 after the simulation of a strong signal decrease applied in the vicinity of an artery. The two inserts show the area of interest before (upper-left) and after (down-right) the application of the signal decrease. (b): Averaged EE over time and over a region of 15 × 15 pixels centered on the simulated heated region, for each tested volunteer. The results obtained with the proposed direct method are compared to those obtained with the minimization method as well as the Horn & Schunck algorithm (Eq. (1)). Numbers above the grey bars recall the number of employed eigenvectors M .

http://bsenneville.free.fr/RealTITracker/

ACKNOWLEDGMENT

The authors thank the European Research Council (project ERC-2010-AdG-20100317, Sound Pharma) and the Dutch Technology Foundation (STW) (project Ontrack) for funding.