Baudouin Denis De Senneville

Charles Mougenot

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

Real time adaptive methods for treatment of mobile organs by MRI controlled High Intensity Focused Ultrasound

Keywords: HIFU, MRI, mobile organ, temperature mapping, motion correction

Real-time adaptive methods for treatment of mobile organs by MRI-controlled high-intensity focused ultrasound

Introduction

High Intensity Focused Ultrasound (HIFU) is a non-invasive technique for local thermal therapy [START_REF] Lele | Production of deep focal lesions by focused ultrasound current status[END_REF][START_REF] Fry | Physical Factors Involved in Ultrasonically Induced Changes in Living Systems: I1 Amplitude Duration Relations and the Effect, of Hydrostatic Pressure for Nerve Tissue[END_REF][START_REF] Stewart | Clinical outcomes of focused ultrasound surgery for the treatment of uterine fibroids[END_REF][START_REF] Furusawa | Magnetic resonance-guided focused ultrasound surgery of breast cancer: reliability and effectiveness[END_REF][START_REF] Pernot | High power transcranial beam steering for ultrasonic brain therapy[END_REF]. It is based on oscillating pressure waves with a wavelength of approximately 1 mm that can be concentrated into a focal point with a size depending on ultrasound frequency, focal length and aperture of the transducer. During the wave propagation across tissue, part of the ultrasound energy is absorbed leading to increase of temperature and eventually to cellular necrosis [START_REF] Dun | Physical mechanism of the action of intense ultrasound on tissue[END_REF][START_REF] Pon | The role of heat in the production of ultrasonic focal lesions[END_REF]. The combination of focused ultrasound and MRI (MR-HIFU) is especially promising since MRI allows a high precision of target localization and it is the only noninvasive imaging technique that has been shown to allow accurate temperature mapping in vivo [START_REF] Cline | MR temperature mapping of focused ultrasound surgery[END_REF][START_REF] Wlodarczyk | Threedimensional monitoring of small temperature changes for therapeutic hyperthermia using MR[END_REF][START_REF] Hynynen | Pre-clinical testing of a phased array ultrasound system for MRI-guided noninvasive surgery of the brain-A primate study[END_REF][START_REF] Kennedy | High intensity focused ultrasound: surgery of the future?[END_REF]. A large range of applications of MR-HIFU has been demonstrated in the clinical field, e.g. in ablation of breast tumors [START_REF] Furusawa | Magnetic resonance-guided focused ultrasound surgery of breast cancer: reliability and effectiveness[END_REF][START_REF] Hynynen | MR imaging-guided focused ultrasound surgery of fibroadenomas in the breast: a feasibility study[END_REF][START_REF] Zippel | The use of MR imaging guided focused ultrasound in breast cancer patients; a preliminary phase one study and review[END_REF] and uterine fibroids [START_REF] Stewart | Clinical outcomes of focused ultrasound surgery for the treatment of uterine fibroids[END_REF][START_REF] Tempany | MR imaging-guided focused ultrasound surgery of uterine leiomyomas: a feasibility study[END_REF][START_REF] Stewart | Clinical outcomes of focused ultrasound surgery for the treatment of uterine fibroids[END_REF] and in preclinical studies on gene therapy [START_REF] Guilhon | Image-guided control of transgene expression based on local hyperthermia[END_REF], and for local drug delivery [START_REF] Mcdannold | Heat-activated liposomal MR contrast agent: initial in vivo results in rabbit liver and kidney[END_REF]. MR-HIFU has shown great potential in improved therapeutic precision, and towards cost reduction.

For optimal treatment, temperature in the targeted volume should be controlled accurately and continuously. To obtain such a control, it has been shown that MR thermometry can be used to provide spatial and temporal temperature feedback control of the HIFU device. Until now, such a control system has been demonstrated for immobilized tissues in vitro and in vivo [START_REF] Vimeux | Real-time control of focused ultrasound heating based on rapid MR thermometry[END_REF][START_REF] Salomir | Hyperthermia by MR-Guided Focused Ultrasound: Accurate Temperature control Based on Fast MRI and a Physical Model of Local Energy Deposition and Heat Conduction[END_REF][START_REF] Quesson | Automatic Control of Hyperthermic Therapy Based on Real-Time Fourier Analysis of MR temperature Maps[END_REF][START_REF] Mougenot | Automatic Spatial and Temporal Temperature Control for MR-guided Focused Ultrasound Using Fast 3D MR Thermometry and Multispiral Trajectory of the Focal Point[END_REF]. However, treatment of mobile tissues (e.g. liver and kidney tumors) would add a major new field of MR-HIFU applications. In order to make this possible, the MR-HIFU system must be adapted to take the motion into account, both with respect to positioning of the focal point (target tracking), as well as to correcting for motion artifacts in temperature mapping. Most tissue displacements during such interventions are due to periodical motions (induced by respiratory or cardiac activity). Such motions may be treated as rigid motion (translation and rotation) but this may be insufficient because of the elasticity of the organ leading to deformation during motion. MR temperature information is obtained using the temperature dependency of the Proton Resonance Frequency (PRF) [START_REF] Hindman | Proton resonance shift of water in the gas and liquid states[END_REF] This equation assumes that the temperature variation T is proportional to the phase variation divided by the magnetic field intensity B 0 , the gyromagnetic ratio , the echo time T E and the frequency temperature constant . Tissue displacements may lead to MR signal phase changes independent of temperature and hence, PRF temperature measurements are highly sensitive to motion artifacts. From the temperature evolution in each voxel the thermal dose can be calculated which is required to predict tissue necrosis [START_REF] Sapareto | Thermal dose determination in cancer therapy[END_REF]. Thus, accurate temperature measurements allow to check the HIFU focal point position and to compute the accumulated thermal dose during the treatment for on-line prediction of tissue ablation. The motion sensitivity currently limits MR-HIFU application to organs with small motion amplitudes or to those that are easily immobilized like uterine fibroid [START_REF] Tempany | MR imaging-guided focused ultrasound surgery of uterine leiomyomas: a feasibility study[END_REF], muscle [START_REF] Hazle | Magnetic resonance imaging-guided focused ultrasound thermal therapy in experimental animal models: correlation of ablation volumes with pathology in rabbit muscle and VX2 tumors[END_REF] and breast [START_REF] Hynynen | MR imaging-guided focused ultrasound surgery of fibroadenomas in the breast: a feasibility study[END_REF].

When intra-scan motion artifacts are limited by using fast imaging techniques and/or synchronization of image acquisition with the periodic motion (e.g. respiratory gating or navigator echo [START_REF] De Zwart | On-Line Correction and Visualiszation of Motion During MRI-Controlled Hyperthermia[END_REF]), image processing can be used to further improve the quality of thermal maps. Some algorithms developed to correct temperature maps have previously been tested in patient studies on the basis of temperature stability in the absence of heating [START_REF] De Senneville | Atlas-Based Motion Correction For On-Line MR Temperature Mapping[END_REF][START_REF] De Senneville | Image processing for on-line reduction of thermometry artifacts[END_REF]. These studies have shown the efficiency of atlas based multiple reference images. The atlas is generated during a preceding "learning phase". Such techniques provide accurate motion estimation and can thus be used to predict the target displacement in periodic motion and to provide information for tracking target tissue. Without such corrections, treatment may be inefficient or may induce unwanted destruction of healthy neighboring tissue. Here, it is shown that a modern, home-built, HIFU system based on a Phased Array ultrasound transducer and rapid electronic displacement of the focal point, in combination with advanced real-time image processing with correction of motion artifacts in temperature mapping, allows tracking of the mobile target during the heating procedure and regional temperature control.

Methods

Focused ultrasound under MR control

A phased array HIFU transducer with a diameter aperture of 96 mm and a focal length of 80 mm, was designed in-house and manufactured by Imasonic SA (Besançon, France). The positioning of the 256 individual elements of 5.8 mm diameter was optimized for minimal R2.17 Real-time PRF based MR-thermometry was performed using RF-spoiled, lipid suppressed (using a 1-1 binomial frequency-selective RF pulse), single-shot gradient echo echo-planar sequences. The following acquisition parameters were used: echo time 40.9 ms, repetition time 1 second, slice thickness 4.5 mm. A single slice with a spatial resolution of 128×128 voxels was reconstructed from 96×96 acquired voxels. With a field of view of 192×192×4.5 mm 3 , the resulting voxel size is 2×2×4.5 mm 3 . Initial data processing was done by the MR acquisition computer. Real and imaginary data were transferred on-line to the HIFU computer in charge of temperature visualization and control of HIFU driver. This HIFU computer is equipped with an Athlon 3.2 GHz XP processor and 1 GB of RAM running under the Windows XP operating system. From the acquired data set an image registration algorithm estimates organ displacement in an area of interest in order to register images and suppress inter-scan motion artifacts in temperature maps. Briefly, an array of reference images results from the initial "learning phase". The procedure makes it possible to define automatically the position of the focal point in the next firing period in order to focus the deposition of energy in the desired tissue area. The R2.9

Evaluation platform

For testing MRI controlled HIFU together with spatial temperature control, a mobile evaluation platform was installed on the patient MRI bed (Figure 1). A plastic bag with fresh pig muscle was installed on the platform. Its displacement was remotely controlled by a mechanical transmission line. The motion period is varied by adjusting the voltage applied to the motor. Also, the motion can be rigid when placed over a sliding support or elastic upon the use of a stopping block. For studying rigid motion, e.g. translation, the position of the transmission line is measured with a millimeter striped bar. Two transceiver optic al diodes detect the position of the striped bar with a precision of 0.5 mm. A microcontroller at 20 MHz (Microchip Technology Incorporated, USA, reference PIC16F688) records the number of graduated displacements with a microsecond temporal accuracy. As soon as a graduation change is detected, timing and position of the target are transferred to the HIFU therapeutic R2.11 R2.8 R2.17

Theory

Since PRF-based temperature maps result from MR signal phase differences, MRI temperature imaging is highly sensitive to intra-scan and inter-scan motion artifacts. Since single-shot acquisition techniques were used in this study, only inter-scan artifacts have to be corrected. In the following, methods are described to correct such artifacts in temperature maps. Then, the position of the focal point is adjusted in order to track the target displacement. Note that data processing has to be performed rapidly to ensure on-line monitoring of temperature evolution. It implies that image processing has to be done within the update time between two successive acquisitions.

Correction of periodical motion artifacts in temperature maps

As temperature is computed from phase difference (see equation [START_REF] Lele | Production of deep focal lesions by focused ultrasound current status[END_REF]), organ motion induces specific inter-scan thermometry artifacts. Indeed, during a therapeutic intervention of a mobile organ, the main parameters involving phase signal modifications are:

1) the temperature variation induced by the heating of the targeted area;

2) the spatial transformation generated by the displacement of the signal source induced by organ motion;

3) the perturbation generated by the modification of the magnetic susceptibility field induced by organ displacements.

It is hardly feasible to separate those three contributions in a single phase image, and therefore the second and third causes should be removed separately. During the intervention, motion field vectors estimated from anatomical images can be used to correct the spatial transformation in phase images (an outline of this processing step is given in Appendix 1). This process improves the precision of on-line temperature measurement by calculating phase differences for pixels corresponding to identical spatial localization. The measured phase variation (and thus the apparent rise of temperature) is also related to the modification of spatial susceptibility effects [START_REF] Salomir | A fast calculation method for magnetic field inhomogeneity due to an arbitrary distribution of bulk susceptibility[END_REF]. The main variations of the magnetic field are located at the In order to correct artifacts generated in temperature maps by periodical organ motion, the perturbation of phase images upon motion is analyzed during a pretreatment procedure.

This step is performed prior to the therapeutic intervention without hyperthermia. An atlas of motion is constructed with MR images acquired during this pretreatment period [START_REF] De Senneville | Atlas-Based Motion Correction For On-Line MR Temperature Mapping[END_REF][START_REF] De Senneville | Image processing for on-line reduction of thermometry artifacts[END_REF]. For that purpose, 50 images are acquired to sample the periodical motion. The reference anatomical image for motion estimation for the image registration algorithm is chosen to be the first one in the temporal series. During the pretreatment step, anatomical images are stored in the atlas with the corresponding phase image. Then, during the thermal therapy, each time a new data set is available, the current anatomical image is compared to anatomical images stored in the atlas using an inter-correlation coefficient. The anatomical image of the atlas with highest similarity is selected, and the corresponding phase image is used as the reference for temperature computation. Finally, the motion field vector estimated from anatomical images is used to correct the spatial transformation in the obtained temperature map.

Adjustment of focal point position taking during periodical organ motion

During an intervention on a mobile organ, the focal point position should be adjusted in order to focus the targeted area initially defined by the radiologist. Target displacements due to periodical motions require a continuous focal point adjustment. Note that the position estimated on the last acquired image is different with respect to the actual tissue position due to a time delay. This delay between the real displacement of the target and the availability of motion information is composed of time periods needed for acquisition, reconstruction, transfer and processing (thermal map computation, motion estimation) of the data. As implementated in this study, the delay is close to 2 s and thus cannot be neglected. If uncorrected, with a typical period of the human respiratory cycle of approximately 5 s, the estimation of the focal point displacement would be close to an opposition phase with respect to the actual target position. Thus, the error in the focal point positioning may increase by a factor of almost two. In addition, the interval of time for data transfer between the MRI and the HIFU computer may be variable resulting in random temporal variations of the estimated period of the motion. For the reasons above, the target displacement must be anticipated with a delay which must be quantified. The criterion used for this anticipation is based on an analysis of the main global motion variation detected in the MR images. This average displacement offers the advantage to be stable against noise because it is computed from all image voxels. A motion atlas is built from motion field vectors computed during a pretreatment step. Each field is associated to the corresponding main global motion.

The motion anticipation requires several steps:

Estimation of the typical period of the motion

An analytic model of the main global motion is defined at the end of the pretreatment step in order to provide during the intervention: 1) a temporal localization of the displacement in the typical period estimated from the current dynamic, 2) a mathematical anticipation of future displacements until the acquisition of the next dynamic.

Due to possible irregularities of the respiratory cycle, an average of the period duration is evaluated during the pretreatment step to define the duration of a typical period. The pretreatment step covers several periods, which are not necessarily identical. To maintain the method compatible with variable respiratory cycle, all sampled periods of the pretreatment step are projected onto an average period. This procedure leads to a densely sampled average period. This main periodical motion M is described with the N th order Fourier decomposition:

( ) ( ) = + = N n n n t n b t n a t M 0 sin ) cos( [3]
In this equation a n and b n correspond to harmonic values. In practice, N=3 is found to be a good compromise to get an efficient typical period model. Higher harmonics correspond to amplitudes lower than 0.1 mm.

Coefficients a n and b n are computed from a least square approximation method with P points acquires during the pretreatment step as follow:

( ) ( ) = = = = = = P i i P i i i n P i i P i i i n t n t n t M b t n t n t M a

Current phase within the estimated motion

To determine the temporal point t i of an image with an estimated main displacement D in the typical period, the equation D=M(t i ) must be solved. However this equation leads to several solutions. A non-determination subsists between increasing and decreasing parts of the cycle. The history of motion is used to solve this ambiguity. For that, an Euclidian distance is computed between the L th last main displacements (in this study L has been set to 5) and the trigonometric polynomial function. This method provides accurate and stable temporal localization of acquired dynamics in the typical period.

Anticipation of the displacement

With the knowledge of time t i and compensating delay measured by the microcontroller, the real target displacement can be deduced from the model M(t). Future displacements are also computed with a temporal resolution of 100 ms. 15 displacements are predicted in order to provide enough anticipated displacements until the treatment of the next dynamic. Displacement is anticipated for a delay of 1.5 s. This provides a margin on transmission jitter of 0.5s taking into account that the dynamic duration is set to 1 s.

With a translational motion along a single axis, knowledge of the main global motion suffices to determine the displacement of the complete target. For complex motions, the image in the atlas with the closest main global motion is used to select the corresponding motion field.

Automatic temperature control in mobile target

Temperature control in a single HIFU focal point

In a single point, the temperature can be adjusted from MR thermal map [START_REF] Vimeux | Real-time control of focused ultrasound heating based on rapid MR thermometry[END_REF][START_REF] Salomir | Hyperthermia by MR-Guided Focused Ultrasound: Accurate Temperature control Based on Fast MRI and a Physical Model of Local Energy Deposition and Heat Conduction[END_REF][START_REF] Quesson | Automatic Control of Hyperthermic Therapy Based on Real-Time Fourier Analysis of MR temperature Maps[END_REF], by a proportional, integral and derivative control using a controller. This technique is based on the differential equation [START_REF] Pernot | High power transcranial beam steering for ultrasonic brain therapy[END_REF] including proportional, derivative and integral terms to minimize the temperature error , The parameter represents the difference between the target temperature and the measured temperature T, the proportional control (second term) is equal to the instantaneous error between the measured and the targeted temperature, the integral control (third term) is defined by the sum of the past temperature errors, and the derivative control (first term) is determined by the present temperature variation in order to reach the future target temperature.The parameter a can be defined by the operator and influences the relative importance of the three contributions in a PID controller. This automatic control ensures a stable convergence to the target temperature.

The proportional, integral and derivative control is also combined with the bioheat transfer equation [START_REF] Dun | Physical mechanism of the action of intense ultrasound on tissue[END_REF] in order to anticipate the tissue reaction. This equation takes into account the thermal diffusion coefficient D, the ultrasound absorption of the tissue and the ultrasound power intensity P.

P T D t T + = 2 [6]
The PID differential equation is respected if the power is set to the value defined by the expression [START_REF] Pon | The role of heat in the production of ultrasonic focal lesions[END_REF] which adjusts correctly the derivative term:

( ) ( ) + + = t C T a T a T D t T P 4 1 2 2 [7]
To control the temperature in a voxel of mobile tissue, the PID feedback loop algorithm is computed from the corrected thermal map which is registered to a central reference position.

With those thermal maps all computations concerning required power for temperature control are carried out in the same way as for static tissue. After required power determination, the position of the focal is adjusted every 100 ms with the anticipation motion algorithm previously described.

Spatial temperature control

To treat a large volume of several milliliters with a small focal point volume of only one cubic millimeter, a proposed solution is to move the focal point along a trajectory [START_REF] Mougenot | Automatic Spatial and Temporal Temperature Control for MR-guided Focused Ultrasound Using Fast 3D MR Thermometry and Multispiral Trajectory of the Focal Point[END_REF][START_REF] Lele | Induction of deep, local hyperthermia by ultrasound and electromagnetic fields: problems and choices[END_REF]. To spatially control the temperature, the amount of applied acoustical power is defined by equation [START_REF] Pon | The role of heat in the production of ultrasonic focal lesions[END_REF] for all points of the targeted volume . However , whatever the type of technology used to displace the focal point (electronic or mechanical), it is technically difficult to emit this defined power simultaneously into each point. To avoid this problem, the R1.9

) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( )
+ + = t r r r r r r F r T a T a T D t t E r r r r r r r 4 2 2 [8]
This energy is delivered successively into every point, choosing the appropriate duration and power [START_REF] Mougenot | Automatic Spatial and Temporal Temperature Control for MR-guided Focused Ultrasound Using Fast 3D MR Thermometry and Multispiral Trajectory of the Focal Point[END_REF]. Points requiring high energy value are heated preferably with a long duration rather than high power intensity for improved safety. Therefore, the feedback cycle time is divided into several durations with a weighting corresponding to the needed energy in order to limit maximal power emission. To optimize the quality of temperature control, the feedback cycle is set to the shortest possible value which is equal to the duration of one dynamic (1 s in this study).

This control method assumes that there is no energy interaction between successive ultrasound beams for the different focal point locations. In a plane perpendicular to the beam axis, this hypothesis is valid when secondary lobes are negligible. Along the beam axis there is an overlap of the acoustic field of the different focal point locations. This overlap tends to produce an elongated heated volume along the beam axis depending on focal point distance and the transducer aperture angle. To avoid having to take into account such overlap in the feedback algorithm, the temperature control has been limited to a coronal slice perpendicular to the vertical beam axis [START_REF] Mougenot | Automatic Spatial and Temporal Temperature Control for MR-guided Focused Ultrasound Using Fast 3D MR Thermometry and Multispiral Trajectory of the Focal Point[END_REF]. In addition, the temperature increase produced in adjacent parallel slices is inferior to the temperature in the controlled region corresponding to a linear (and not a surface) heating trajectory.

The spatial feedback temperature control is also computed with corrected and registered thermal maps. The trajectory provided by this feedback algorithm is modified in order to force focal points to track the target. Each point of this trajectory is split in a succession of several points with durations close to 100 ms. After splitting, those points are translated with a shift defined by the anticipated motion. As a result, the energy defined by the spatial temperature control algorithm is deposited with a high sampling rate.

Results

To test the efficiency of the proposed method, ex-vivo muscle was heated with a HIFU transducer placed below the evaluation platform during successive applications of periodical rigid and elastic motion.

Periodical rigid motion

Ex-vivo muscle was displaced with a periodical translational motion of 14 mm with a duration of 5.6 s. The muscle was heated with HIFU emitting 24 Watt acoustical, during 1 minute. Figure 2 gives a comparison of temperature maps (during the heating) and the thermal dose maps (at the end of the experience) between studies realized with and without thermal artifacts and focal point position correction. Temperature errors are reduced from 40 °C to less than 1 °C when the thermal correction is applied making focal point localization and thermal dose computation possible. Without focal point position correction, the heating area is spread along 14 mm length and the induced heating does not lead to tissue ablation as temperature increase is insufficient. On the other hand, with correct focal positioning, the heating pattern is circular and the energy is deposited over the same tissue area leading to necrosis.

Periodical elastic motion

The more frequent and complex case to correct for is periodical elastic motions. To mimic those displacements, ex-vivo muscle tissue was submitted to elastic periodical motion with a period duration close to 5.6 s. The muscle has been heated by HIFU emitting a power of 100 electrical watts during 1 minute. Figure 3 shows anatomical images obtained at different positions in the cycle with the corresponding field vector displacements estimated with a centered reference image. Figure 4 The heating procedure with the focal point position correction, described previously, provides a temperature increase with a constant power deposition in the same target point despite tissue motion. Those correction methods also provide good temperature maps with a quality close to the maps acquired without motion. Thus, the corrected temperature maps may be used to control temperature. PID feedback loop algorithms were applied to an ex vivo muscle tissue undergoing a periodical translation.

Figure 5 shows a comparison between two heating procedures realized with spatial temperature control of immobile tissue (5a and 5c) or undergoing periodical rigid motion (5b and 5d). The amplitude of the motion has been set to 8 mm along the Z axis in order to maintain the possibility to move the focal point 9 mm along the perpendicular axis X.

According to those amplitudes, the focal point does not deviate more than 6 mm from the central point. With a voxel width of 1.5 mm, the temperature control was achieved over 7 central voxels to obtain a linear heating segment of 9 mm. Light grey, dark grey and black solid curves on graphs 5a and 5b represent minimal, average and maximal temperature, respectively, for those 7 voxels compared to dashed curves representing the targeted temperature. Light grey, dark grey and black solid curves on graphs 5c and 5d represent temperature measures along this segment at three different time points during the temperature increase. Light grey, dark grey and black dot curves describe the corresponding target temperature increase.

The actual temperature increase corresponds to the targeted temperature increase with an accuracy of 0.5 °C over the entire heating segment. The average difference between the maximal and the minimal temperature was 1.2 °C for the study with immobile tissue and 1.3 °C with mobile tissue. With the motion tracking technique, tissue displacements did neither induce any significant noise in temperature maps nor lead to any measurable loss of precision in temperature control. Figure 6 shows a comparison of temperature maps acquired at 141 s and 215 s obtained ex-vivo with immobile and muscle undergoing periodical rigid motion in the direction top-to-bottom of Figure 6. The elongated segment of 9 mm along the X axis appears clearly as if the tissue moved along the perpendicular direction with an amplitude of 8 mm. The straight segment shapes confirm that focal point motion tracking was achieved correctly.

Since the focal point position is adjusted as a function of the tissue displacement, the temperature control accuracy only depends on feedback algorithm performance. Because the differential equation [START_REF] Pernot | High power transcranial beam steering for ultrasonic brain therapy[END_REF] of the PID controller gives a dual solution of a/2, the response time of R1.10 the system is equal to 2/a. To avoid temperature control instability induced by system latency around 2 s (measurement, transfer and analysis delay), the response has been set to the duration of 8 s.

Precision of the tracking method proposed for periodical motions

To evaluate the performance of the correction method, it is interesting to compare the estimated displacement from MR images and the anticipated displacement with the real displacement of the target measured by a microcontroller via a striped bar placed on the transmission line.

Approach used for rigid periodical motions correction

For stable motion compensation, it is necessary to update the anticipated displacement as often as possible. The anticipated displacement is updated at each acquired dynamic with a measurement of the delay to compensate. This delay is mainly composed of the dynamic acquisition duration (1 s), the transmission and the processing delay which varies from 0.3 s to 1.1 s, with an average value of 0.9 s. This variation of the availability of motion information depends on the operating system used. In the example shown this was carried out under Windows XP (not a real time Operating System). The delay was measured directly with a microcontroller keeping time between the beginning of the acquisition of the last dynamic (external pulse provided by the image reconstructor) and the end of the corresponding processing. Figure 7 describes results obtained with an anticipated motion, updated at each dynamic based on this measured delay. Those results correspond to the study with a periodical translational motion reported in Figure 2c and2f. The method is shown to be efficient since the standard deviation between the measured and the anticipated motion remains constant close to 0.3 mm in the experiment reported in Figure 2. The anticipated motion is a continuous curve except for infinitesimal discontinuities matching the periodical target motion.

To evaluate precisely the accuracy and the reproducibility of the target position anticipation process, this experiment was repeated 14 times. The average error and standard deviation over all studies was 0.33 ± 0.03 mm, respectively. The minimal and maximal errors evaluated over one period were 0.26 mm and 0.45 mm, respectively.

Approach used for correction of elastic periodical motions

For elastic deformations, the main displacement was insufficient to characterize the estimated motion of the target. Moreover, the displacement measured with the transmission line is not enough to characterize the target point in the tissue and can not be compared to the motion tracking for elastic motions reported in Figure 3. To correct for elastic motions, the main displacement was used as a criterion to select appropriate field vectors in the atlas.

Therefore, to quantify this lack of performance, the anticipated motion obtained when the atlas technique was applied to a periodical translational motion similar to the previous one was compared to the motion measured by the striped bar. The anticipated displacement was sampled with 50 values corresponding to the 50 field vectors stored in the atlas. This sampling of the anticipated motion decreased slightly the accuracy of the motion determination. The error between the real and the anticipated motion was close to 0.4 mm (average over 4 experiments) as compared to 0.3 mm obtained previously using rigid body motions.

Discussion

The use of phased array transducer technology allows very fast electronical displacement of focal point. The transducer with 256 elements of 6 mm diameter used in this study allowed a displacement of 15 mm at 1.5MHz. This amplitude offers the possibility to correct most of the displacements in organs under respiratory motion. The described motion correction methods can be directly generalized for three-dimensional motions. Results described in this study are two-dimensional because of limitations imposed by transmission of data between MRI computer and HIFU computer.

Signal localization errors may be caused by instrumental imperfections and/or limitations. Some localization errors can also be induced by the imaged object, like susceptibility artifact for paramagnetic material or chemical shift for a specific water and fat composition. As a result, geometric distortions in modulus and phase images modify anatomical images and estimated motions. Those artifacts can be limited with an optimal MR acquisition sequence or corrected with specific correction processes [START_REF] Schmithorst | Simultaneous correction of ghost and geometric distortion artifacts in EPI using a multiecho reference scan[END_REF].

Focal point position correction upon organ displacement is necessary for treatment of a mobile target. However, in case of MRI guidance a delay in execution of the correction F o r P e e r R e v i e w cannot be avoided. Before an image is available, there is an acquisition, transmission and calculation delay close to 2 s. During this duration the position of the focal point is not yet updated. For periodical motions, the correction technique described is limited to motion with a period duration larger than the acquisition duration. According to the Shan non-Nyquist theorem, the sampling requires at least 2 dynamics per period to reconstruct a typical period model . In the studies shown, the motion period duration was 5.6 s and the dynamic acquisition duration was 1s, sufficient to obtain a good period reconstruction.

Displacements applied to correct motion artifacts are limited to displacements detected during pretreatment step. This limitation of motion sampling may be avoided using interpolation of detected displacement.

Analysis of the precision of the focal point repositioning

Without displacement correction, the average error of the focal point position was measured experimentally at 4.8 mm. A theoretical estimation of this error can be calculated with a first harmonic approximation of the real motion. Equation [START_REF] Stewart | Clinical outcomes of focused ultrasound surgery for the treatment of uterine fibroids[END_REF]

= = = + = [9]
with c 1 indicating the amplitude of the first harmonic and using the definition of a and b as in Equation 3. The theoretical error induced by this motion is:

( ) ( ) mm c dt t w c T T 6 , 4 2 cos 1 1 0 2 1 correction without = = = [10]
When the focal point was repositioned with the displacement estimated from the last available image but without taking into account the anticipated displacement during the delay following the last available image, the average experimental error was 7.54 mm. The theoretical estimation of this error might be quantified by supposing that the delay to be compensated is constant (d C =1.9 s). The difference in phase induces is equal to C :

°= = = 120 1 , 2 ' 2 rad T d C C [11]
The theoretical error induces by this phase difference is: When the delay to be compensated is correctly measured, the experimental error measured is 0.3 mm with rigid motion algorithm and 0.4 mm with elastic motion algorithm. Theoretically, there is no phase difference between the anticipated and the real displacement, therefore this error should be 0. Table 1 shows a comparison of the theoretical and experimental standard deviations obtained with each approach for the periodical translation motion. The anticipated motion is shown to provide an accurate focal point repositioning, 14 times better as compared to positioning without motion correction.

( ) ( ) ( ) 8 2 sin 2 cos cos 1 1 0 2 1 1 mm c dt t w c t w c T T C n nticipatio without a = = + = [12]

Stability of the displacement correction in case of motion periodicity variations

The duration between the localization of the current dynamic in the typical period model and the required anticipation time is composed of:

-The localization process of a dynamic within the motion cycle. This step required a windowing of L dynamics. The average delay induced is equal to (L-1)/2 dynamics duration.

-The delay to compensate which is composed of the acquisition, transmission and processing duration of one dynamic.

-The anticipation of the trajectory until the next dynamic acquisition. The average of this delay is equal to half a dynamic duration.

This duration corresponding to the system response delay d R is thus equal to 2+1.9+0.5=4.4 s.

The periodical motion hypothesis is not perfectly respected in practice because of variations in physiological cycle. Therefore, the period duration may change during the intervention.

This induces an error in motion anticipation. If the period motion duration is equal to T, the anticipated motion is calculated with a phase shift equal to = 2 •d R /T. On the other hand, if the period duration becomes T', the phase shift required is ' = 2 •d R /T'. However, the anticipated motion is still computed with a phase shift instead of '. Therefore, the phase difference between the anticipated motion with a period duration T and the real motion with a period duration T' is :

! = = T T d R 1 ' 1 2 ' [13]
In this case, the standard deviation between the real and the anticipated motion is: To study the robustness of the proposed correction method, several tests were carried out with a period duration motion change of ±1 s with respect to the initial period of 5.6 s.

Figure 8 compares the experimental standard deviation measured with a modification of the period duration and the theoretical standard deviation as described in Equation [START_REF] Tempany | MR imaging-guided focused ultrasound surgery of uterine leiomyomas: a feasibility study[END_REF]. It can be seen that a standard deviation of 1.5 mm (corresponding to the width of a voxel) is obtained when the period variation does not exceed ±7%. The stability of this method against variation of period duration is proportional with the following parameters: i) the system response delay, ii) the motion amplitude, iii) the estimated period during the pretreatment step. In the general case of periodical motion with n harmonics the accuracy of the anticipated displacement is described by this expression where c n indicates the amplitude of the nth harmonic, and the definition of a and b as in Equation 3.

Precision of reported results

The uncertainty of reported results depends on the estimated motion method used and the accuracy of the evaluation platform. The two main sources of uncertainties induced by the platform are:

-The striped bar details the position of the transmission line with an accuracy of 0.5 mm. However, intermediate values are computed with a linear interpolation from the detected position.

-The mechanical connection between the transmission line and the ex vivo muscle is not perfectly rigid because of external friction. Those frictions may cause an estimated error of 0.5 mm of at each change of motion direction.

In fact, the actual standard deviation between the anticipated and the real motion is expected to be significantly lower than the value of 0.3 mm measured with this evaluation platform.

Conclusion

Correction methods for displacement artifacts in temperature maps of mobile targets may reduce temperature errors from several hundred degrees to values lower than one degree. The techniques proposed here provide correction of focal point position for periodical complex motions (rigid or elastic). The correction of focal point position is accurate with a standard deviation error lower than 0.3 mm for the studies reported here. The temperature increase obtained for MR-HIFU with this correction is more efficient and correctly localized despite target displacements. Also, it is shown here that MR-HIFU heating with motion tracking is possible with spatial temperature control. This may provide tissue ablation without damaging neighboring tissue during motion.

The presented methods have been tested for two dimensional motions and can be generalized for three dimensional motions. The proposed techniques appear stable enough to treat organs like kidney or liver with MR-HIFU for a patient under respiratory control. For patient under free breathing, the system response delay should preferably be further decreased.

Real time information like respiratory gating, navigator echo or ultrasound echo [START_REF] Pernot | 3D Real-time Motion Correction in High Intensity Focused Ultrasound Therapy[END_REF] may be used instead of main global motion as a selection criterion for displacement field in the atlas. With further acceleration of motion tracking data, the stability is expected to be improved further because the anticipation delay will be reduced correspondingly. [START_REF] Friston | Spatial registration and normalisation of images[END_REF]. In a second step, a hierarchical approach of Horn and Schunck algorithm (37) provides a good estimation of local organ displacements because of the regularity constraint assuming that motion field vectors have similar values for adjacent pixels matching real organ motion. Note that displacement estimation based on MR images allows sub-pixel accuracy as the registration process is computed with a large number of pixels.

Such image registration required 200 ms of computation time on our HIFU computer for one image of resolution 128×128 pixels. Further details can be found in references [START_REF] Maintz | A survey of medical image registration[END_REF][START_REF] Barron | Performance of Optical Flow Techniques[END_REF].

  according to the following equation for phase changes in RF spoiled, gradient echo sequences (23,24):

  lobes and compactness of the transducer[START_REF] Gavrilov | A theoretical assessment of the relative performance of spherical phased arrays for ultrasound surgery[END_REF]. The system is able to focus an ultrasonic wave at 1.5 MHz in a focal point with a size of 0.8×0.8×3.5 mm³ (-3 dB level) according to acoustic field simulation based on Rayleigh integration. The acoustic pressure and the position of the focal point are adjusted by choosing the amplitude and the phase of electrical signals. The position of the focal point is adjusted electronically over a volume of 15×15×30 mm³ at a rate of 100 ms. The electrical to acoustic signal conversion efficiency of the transducer is 37 % with an impedance adaptation to 50 Ohm. In addition, the loss of the transmission line amounts to 26 %. This efficiency has been measured for each channel with radiation force pressure measurement using a Sartorius balance (precision of 0.1 mg, model CP324S, Goettingen, Germany) and an absorbing target provided by Imasonic (Besançon, France) with an elastic rubber material filled with metallic powder in order to fully absorb the acoustic energy.Experiments were performed on a Philips 1.5 Tesla Intera clinical scanner (Philips Medical Systems, Best, the Netherlands). A sample holder allowed the US beam to propagate through a water-filled volume from transducer to sample. A receive-only surface coil (92 mm diameter) was positioned parallel to the transducer plane. Demineralized and degassed water has been used at the water-sample interface to avoid air bubbles. The oxygen gas content was 4 mg/l (ppm) as measured with an oxymeter (Consort, France, reference C535).

  derived from previously described procedures and is further elaborated below. Periodical motions also require accurate timing measurements in order to anticipate the target position due to periodic tissue displacements. A microcontroller at 20 MHz (Microchip Technology Incorporated, reference PIC16F688) keeps track of the start of the MR pulse sequence with a microsecond temporal resolution. This information is transferred to the HIFU therapeutic computer in a few milliseconds. The position of HIFU focal point and desired power are transferred via optical fibers to an in-house designed 256 channel generator/amplifier, built by Radio Concept Développement et Communication (St Martin de Hinx, France). This generator produces sinusoidal waves using direct digital synthesizer provided by Analog Device (Northwood, England, reference AD9835). Each electrical sinusoidal is amplified with a dual differential linear amplifier manufactured by Texas Instruments (Dallas, United States, reference THS6012). Output signals are isolated with a transformer and filtered using a low pass filter and a 64 MHz reject filter. The electrical signals are then transmitted to the ultrasound phased array transducer via an impedance matching box located next to the magnet cryostat. A specific phase for each channel is produced by direct digital synthesizer components (accuracy 1°) in order to focus the beam at the desired location. The system is designed for a maximum of 180 Watts acoustical power. At the surface of the transducer composed of 256 elements of 26.4 mm², the maximal available acoustical power is 2.7 W/cm² .

  RS232 connection in few milliseconds. Those real-time target position data are used to quantify the accuracy of the estimated and anticipated displacements derived from MR images.

  with different susceptibility. The influence of susceptibility variation on the measured phase signal depends on imaging parameters (spatial homogeneity of the static magnetic field in the MRI, echo time…). Modeling of inhomogeneous magnetic susceptibility field variations upon motion in-vivo is difficult and its feasibility remains unproven for clinical application. Approaches described in this study allow to correct displacement-related errors in PRF based MR thermometry for periodic motions avoiding explicit modeling of the susceptibility field.

  energy E deposed during one feedback cycle t F into the point r is defined by equation[START_REF] Cline | MR temperature mapping of focused ultrasound surgery[END_REF].

(

  provides a comparison between temperature maps (during the heating) and thermal dose maps (at the end of the experience) obtained with and without thermal and focal point position correction.Temperature errors without thermal correction can reach 300 °C in one minute. Phase variations greater than 2 between two dynamics induce an accumulated thermal error which increases during the heating. Similar to experiments described above, thermal dose maps obtained without thermal correction are unusable. Without focal point position correction, the heating area is spread along a 12 mm trajectory and the induced heating does not lead to tissue ablation as temperature increase is insufficient. With correct focal positioning, the heating pattern is circular and the energy is deposited over the same tissue area leading to necrosis.Temperature control with motion tracking

Figure 1 .

 1 Figure 1. Mechanism used to displace the target in a cone-shaped piece of meat. The HIFU transducer is placed in the center of the magnet, and the motor at the end of the patient bed.

Figure 2 .

 2 Figure 2. Temperature variation maps (a, b, c) and thermal dose maps (d, e, f) obtained exvivo during heating of pig muscle undergoing a periodical translation: a and d: without thermal and focal point position correction, b and e: with thermal correction but without focal point position correction, c and f: with thermal and focal point position correction.

Figure 3 .

 3 Figure 3. Anatomical images (a, b, c, d, e) acquired at different time periods during the periodic elastic motion and corresponding displacement field vectors (f, g, h, i, j) required for registering the anatomical image at the reference position.

Figure 4 .

 4 Figure 4. Temperature variation maps (a, b, c) and thermal dose maps (d, e, f) obtained exvivo during heating of pig muscle undergoing an elastic periodical motion: a and d: without thermal and focal point position correction, b and e: with thermal correction but without focal point position correction, c and f: with thermal and focal point position correction.

Figure 5 .

 5 Figure 5. Spatial temperature control along a 9 mm line on immobile tissue (a, c) or tissue undergoing rigid periodical motion (b, d): a and b: temporal temperature measurements, c and d: spatial temperature measurements.

Figure 6 .

 6 Figure 6. Temperature variation maps at 141 s (a, b) and 215 s (c, d) obtained ex-vivo during heating of pig muscle: a and c: with immobilized muscle, b and d: muscle tissue undergoing periodical rigid motion.

27 Figure 7 .

 277 Figure 7. Real (dotted line), estimated (dashed line) and anticipated displacement (solid line) measured with a periodical translational motion during 140s (a) or 20s after 2 minutes (b).

Figure 8 .

 8 Figure 8. Stability of the focal point position correction method as a function of the variation of the duration of the periodical motion for an anticipated delay of 4.4 s.

Figure 1 .Figure 2 .Figure 3 .Figure 4 .Figure 5 .Figure 6 .Figure 7 .Figure 8 .

 12345678 Figure 1. Mechanism used to displace the target in a cone-shaped piece of meat. The HIFU transducer is placed in the center of the magnet, and the motor at the end of the patient bed.

Table 1 .

 1 Error of focal point position for different types of correction in mm.

	R1.11
	R1.12

Page 25 of 37 Magnetic Resonance in Medicine Magnetic Resonance in Medicine

	List of figures
	F o r
	P
	e
	e r
	R e v
	i e
	w

Magnetic Resonance in Medicine Magnetic Resonance in Medicine

Appendix 1: Estimation of tissue displacement during motion

The proposed approach consists of using image processing techniques to estimate online organ displacements from anatomical images [START_REF] De Senneville | Atlas-Based Motion Correction For On-Line MR Temperature Mapping[END_REF][START_REF] De Senneville | Image processing for on-line reduction of thermometry artifacts[END_REF]. Image registration allows organ displacement estimation with an excellent spatial resolution as compared to other existing approaches [START_REF] Pernot | 3D Real-time Motion Correction in High Intensity Focused Ultrasound Therapy[END_REF]. The objective is to relate the coordinate of each part of tissue in the image to register with the corresponding one in a reference image. Here, the general principles are summarized.

In general, 3D vector field maps are preferred for this purpose. However, in practice it is difficult to acquire on-line 3D isotropic images because of the technical limitations and SNR associated with the fast 3D acquisition sequences. The approach used consists of estimating motion on 2D images generated by organs moving in 3D space. Image sections are oriented to include the principal axis of the organ displacement so that apparent movement is the best approximation of the real displacement. However, the presented algorithms may be implemented in the 3D case in order to estimate motion field between two volumes when technological progress will allow sufficient spatial and temporal resolution. (  )

where M is a criterion determining the accuracy of the registration.

Registration efficiency is related to the nature and domain of the transformation to detect (global [START_REF] Friston | Spatial registration and normalisation of images[END_REF][START_REF] Wolberg | Robust Image Registration Using Log-Polar Transform[END_REF] or local [START_REF] Schunck | Determining optical flow[END_REF][START_REF] Lucas | An iterative image registration technique with an application to stereo vision[END_REF] transformation). A compromise has thus to be found between the permissiveness of the spatial transformation T and the robustness of the estimation process. A well-adapted technique consists of using results obtained with a global transformation estimation as a starting point for a local transformation estimation. In this study, a global affine transformation is estimated in a first step, using a differential approach R2.14

Magnetic Resonance in Medicine Magnetic Resonance in Medicine


Karsten Østergaard Noe

Baudouin Denis De Senneville

Vindelev Ulrik

Kari Tanderup

Thomas Sangild Sørensen

Vindelev Elstrøm

Acceleration and validation of optical flow based deformable registration for image-guided radiotherapy

Introduction

Organ deformation during successive image guided radiotherapy (IGRT) sessions represents a significant challenge to optimal planning and delivery of radiation doses. To facilitate a more precise conformation of doses to the tumor and sparing normal tissue, multiple 3D cone beam CT (CBCT), conventional CT or MRI datasets can be acquired during the treatment course. To update existing dose plans based on this newly obtained data, the data must be correlated to a reference dataset. For this task a registration method is needed. A per pixel-based (deformable) registration technique is required to fully account for the non-homogeneous deformation throughout the volume. However, such methods are much more time consuming to perform compared to methods estimating homogeneous (rigid or affine) transforms. Moreover, if the deformable registration is not regularized sufficiently it can result in physically non-plausible deformations causing significant errors in subsequent dose calculations. In this paper we evaluate the accuracy of two fully automated image based deformable registration methods driven by the concept of optical flow [START_REF] Horn | Determining optical flow[END_REF], [START_REF] Cornelius | Adapting optical-flow to measure object motion in reflectance and X-ray image sequences[END_REF]. These methods were chosen based on our previous experience with registration of 2D MRI for online MR temperature monitoring [START_REF] De Senneville | An optimised multi-baseline approach for on-line MR-temperature monitoring on commodity graphics hardware[END_REF]. As a preliminary validation study we report the result of using one of these methods for registering 4DCT lung acquisitions and the other for registering head and neck CBCT and conventional CT acquisitions.

To compute such complex registrations in a clinically acceptable time frame, we implemented the two algorithms in parallel on a commodity graphics processing unit (GPU), an emerging platform for general purpose computation. The methods used in the study presented here are easily parallelizable making them ideal for GPU implementation. Furthermore, the running times of these methods are relatively short even in a non-accelerated version meaning that a significant acceleration will allow us to do deformable registration in very short time frames. Several other deformable registration methods have already been demonstrated with significant speedups on GPUs [START_REF] Strzodka | Image registration by a regularized gradient flow -a streaming implementation in dx9 graphics hardware[END_REF], [START_REF] Noe | GPU accelerated Viscous-fluid deformable registration for radiotherapy 16th Medicine Meets Virtual Reality Conference[END_REF], [START_REF] Sharp | Gpu-based streaming architectures for fast cone-beam ct image reconstruction and demons deformable registration[END_REF]. Consequently, this paper also includes an evaluation of the running times of the two algorithms.

Materials and methods

Optical flow based registration

The process of estimating optical flow means finding a quantitative measure of how image intensity information has changed between two images. Technically both images are regarded as part of one mathematical function where spatial changes have occurred in the time between acquisitions transforming one image into the other. The optical flow is a vector field consisting of the changes in space coordinates. These vectors can be thought of as 'optical velocity' vectors showing the direction of image intensity flow. We focus on two optical flow based methods for deformable image registration; a 3D version of the Horn & Schunck algorithm [START_REF] Horn | Determining optical flow[END_REF] and the extension of this algorithm by Cornelius & Kanade [START_REF] Cornelius | Adapting optical-flow to measure object motion in reflectance and X-ray image sequences[END_REF] to handle intensity changes that do not arise as a direct consequence of geometric motion, i.e. intensity variation due to physical properties of the acquisitions themselves.

The Horn & Schunck method is based on an assumption of preserved image intensity in the two 3D images to be registered. This means that it only works for registering images of the same modality and only for images with consistent grey values when multiple image sets are compared. It is also assumed that the deformation is smooth. This is in general a valid assumption for soft tissue deformation. The Horn & Schunck method is thus suitable for registering successive CT images due to the reproducibility of Hounsfield Units for this modality.

When registering MR to MR or CBCT to CT the assumption of intensity preservation is no longer valid. In the CBCT modality the Hounsfield Units are affected by the larger contribution from x-ray scatter. Also the design of the detector and the image reconstruction algorithm used has an impact on Hounsfield unit reproducibility for CBCT. To facilitate handling of intensity differences, Cornelius & Kanade extended the original algorithm thus enabling it to tolerate some deviation from the assumption. In their work it was further assumed that the non-motion-related intensity differences are smoothly varying in space.

Both our implementations utilize a multi-resolution approach. This means that the organ deformation is first approximated on low resolution versions of the 3D acquisitions to be registered. The result of this coarse registration is then used as a starting point for a registration at a higher resolution. This continues until the deformation has been approximated at the highest resolution. This strategy enables us to systematically handle modes of deformation at different scales.

GPU based registration

The very specialized parallel hardware architecture of modern GPUs enables them to perform a vast number of numerical calculations in a short time frame allowing them to outperform the CPU in normal computers for a number of applications. The degree of GPU acceleration of an algorithm attainable depends on how suitable the algorithm is for being split into a large number of small computations that can be run simultaneously. The registration methods we focus on here are very suited for such a parallelization, and thus we have utilized the computational capabilities of a GPU for accelerating the computations required in the Horn & Schunck and Cornelius & Kanade methods. The GPU implementations are based on the CUDA framework from Nvidia [START_REF]NVIDIA CUDA: Compute unified device architecture -Programming guide[END_REF].

Image material

The Horn & Schunck registration has been evaluated on the POPI-model which is a 4D thorax virtual phantom [START_REF] Vandemeulebroucke | The POPI-model, a point-validated pixel-based breathing thorax model[END_REF]. It consists of 10 CT data sets of resolution 482×360×141 which have been acquired at different breathing phases during a single breathing cycle.

The images were acquired at the Léon Bérard Cancer Center, Lyon, France. In each data set corresponding to a breathing phase 41 landmark points have been manually identified, and these points are used for our validation. The voxel spacing of the acquisitions is 0.98×0.98×2.0 mm 3 .

The image material used for validating the Cornelius & Kanade method is a series of 6 CBCT images of a head and neck cancer patient and a conventional planning CT image acquired at the Department of Oncology, Aarhus University Hospital. The scans have been conducted weekly during the treatment course starting at the first fraction. The CBCT images are of dimensions 512×512×51 with a voxel spacing of 0.47×0.47×3.0 mm 3 , while the conventional CT image is of dimensions 512×512×55 with a voxel spacing of 0.78×0.78×3.0 mm 3 . Validation of bone alignment is based on 6 landmark points in each 3D data set. These points have been manually positioned prior to registration at easily identifiable positions on the bony anatomy of the cervical vertebrae and the base of skull representing clinically relevant match points. Positioning of these points is illustrated in figure 1.

Registration validation studies

Three series of registration experiments have been carried out:

1) CT to CT registration using the Horn & Schunck registration method:

Following the convention from the POPI initiative all images from the 4D data set (at time phases numbered 0, 2, 3, 4, 5, 6, 7, 8 and 9) have been registered to the reference image at phase 1.

2) CBCT to CBCT registration using the Cornelius & Kanade registration method:

The CBCT images numbered 2 to 6 have been registered to CBCT image 1.

3) CBCT to CT registration using the Cornelius & Kanade registration method:

The CBCT images have been registered to the planning CT image.

For the Horn & Schunck method a reference CPU based implementation has also been implemented allowing us to compare the computation times between the CPU and GPU versions.

Image preprocessing

The CBCT images were processed using GREYCstoration image denoising filter [START_REF] Tschumperlé | Fast Anisotropic Smoothing of Multi-Valued Images using Curvature-Preserving PDE's[END_REF]. As the publicly available implementation of this filter works in two dimensions, an in-house program has been used that simply filters each image slice independently. The effect of the filter is to remove noise (and in some cases artifacts from the CBCT reconstruction) while preserving the edge contrast between different kinds of morphology. See figure 2 for an example.

The registration methods we present in this paper are designed to estimate the detailed deformation of morphology. If there is global displacement of patient position (that is translation and/or rotation) between two images it is necessary to do a rigid alignment of the images before the deformable registration in order to supply the method with a suitable starting point for estimation of organ deformation. The rigid registration method we use is based on the Insight Registration and Segmentation Toolkit (ITK). The measure used to compare images is based on mutual information. Input images are filtered using a threshold filter so that only bone morphology is included in the rigid registration. Again a multi-resolution approach is taken.

For validation study 1 no image preprocessing was required as the images were already rigidly aligned. In validation studies 2 and 3, a rigid registration was required. A bounding box corresponding to the physical extent of CBCT image 1 has been cut out of the planning CT image and resampled to the same resolution as the CBCT images.

Results

Validation study 1: Registration of the POPI 4DCT data set

The registration accuracy, evaluated on the target registration error (TRE) of landmark positions, is summarized in table 1. The distances are calculated as the Euclidian length of 3D vectors. Original average landmark distance was 3.5 mm ± 2.0 mm. After registration, this average distance was equal to 1.1 mm ± 0.6 mm.

A visualization of the registration result can be seen in figure 3. In this visualization the source image is shown in a reddish color while the reference image is shown in a bluish color. Where the images align a gray scale image emerges. In the unregistered case on the left blue and red areas can clearly be seen indicating that the morphology is not aligned. In the registered case to the right these colored areas have almost disappeared indicating that the images have been successfully registered.

Validation study 2: CBCT to CBCT registration

In table 2 the alignment error of the landmark points positioned on the bony anatomy is summarized both before registration, after the rigid registration, and after the deformable registration. Original average landmark distance was 5.8 mm ± 1.1 mm. After the rigid registration, this average distance was equal to 1.8 mm ± 1.0 mm and after the deformable registration it was 1.6 mm ± 0.8 mm.

The result of registering CBCT image 3 to CBCT image 1 is visualized in figure 4. The images depicting the results of the rigid registration show an acceptable alignment of most bony anatomy, but it can be seen that the soft tissue and the area surrounding the oral cavity is not aligned. Improved alignment is obtained as a result of the deformable image registration as shown in the visualization. In figure 5 a visualization of the computed transformation can be seen showing that the deformation is smooth.

Validation study 3: CBCT to planning CT registration

The average alignment errors after the CBCT to CT registrations are found in table 3.

After the rigid registration, this average distance was equal to 2.2 mm ± 0.6 mm and after the deformable registration it was 1.8 mm ± 0.6 mm.

The result of registering CBCT image 6 to the planning CT image can be seen in figure 6. Again it can be seen that even though the skull and spine is aligned after the rigid registration, deformable registration is needed to account for changes in jaw positioning and deformation of soft tissue.

Time consumption

On an Intel Core 2 CPU at 2.4 GHz the Horn & Schunck registration used for the each 3D image in the POPI dataset in validation study 1 takes 30 minutes. On an Nvidia Geforce 8800GTX GPU in the same machine each registration takes 37 seconds, making the GPU version 48.7 times faster.

For the 3D Cornelius & Kanade method we did not write a CPU reference implementation. Subsequently we do not know the exact difference in processing time between CPU and GPU. However we expect the acceleration of this method to be somewhat smaller than for the Horn & Schunck method due to a less efficient use of socalled shared memory on the GPU. Each 3D registration of the CBCT images in studies 2 and 3 using the Cornelius & Kanade method takes 64 seconds.

Discussion

When specifying landmark points the limiting factor on the accuracy is the slice thickness. This is because a point may be between two slices which makes it hard to identify. In the light of this we consider the registration accuracy results in validation study 1 very acceptable as the mean landmark error is well below the slice thickness. This accuracy is comparable to results for the Demons algorithm previously reported in the POPI initiative [START_REF] Tschumperlé | Fast Anisotropic Smoothing of Multi-Valued Images using Curvature-Preserving PDE's[END_REF].

In validation study 2 and 3 we used landmark points to track the registration of clinically relevant points on the cervical vertebrae and the base of skull. For this preliminary study we did not do a dedicated evaluation of the error in landmark point identification but we estimate that the error in each landmark point position may be as high as 2.0 mm. The magnitude of this error is mainly due to the slice thickness of 3 mm.

It is interesting to note that although the rigid registration was done on bone morphology the mean landmark error on bone morphology is reduced by the deformable registration.

In all registrations the mean error is smaller than the slice thickness. As demonstrated in figures 4 and 6 a rigid registration is not sufficient in describing the geometrical difference between the images. However these geometrical differences have been substantially reduced by the deformable registration.

Based on these studies we are optimistic that the Cornelius & Kanade method is suitable for registering head and neck CBCT images from a series of radiotherapy treatments to the planning image. Hopefully this will allow us to compensate for the unreliable Hounsfield units by using the inverse transformation of the one found in validation study 3 to map the Hounsfield units from the planning CT to each CBCT. This will make it possible to evaluate the doses delivered in the treatment fraction corresponding to each CBCT acquisition by doing a dose calculation on the corrected CBCT image. This has previously been suggested by Yang et al. [START_REF] Yang | Evaluation of on-board kV cone beam CT (CBCT)-based dose calculation[END_REF]. When all CBCT images are registered to the same geometrical reference system, it will then be possible to evaluate the actual accumulated dose from a series of fractions for comparison with the planned doses. This can be done by deforming the dose distributions from each CBCT to the geometrical reference frame constituted by the planning CT [START_REF] Daniel | Parotid Gland Dose in Intensity-Modulated Radiotherapy for Head and Neck Cancer: Is What You Plan What You Get?[END_REF]. Furthermore an accurate registration makes it possible to do automatic segmentation by transferring segmentations from the planning CT. Currently a larger head and neck CBCT registration study is being planned in which landmarks in bone as well as soft tissue will be used for accuracy evaluation. Also a smaller slice thickness will be used allowing us to more accurately positioning landmarks.

In the pursuit of online IGRT, the performance of the required image processing in a sufficiently short time frame constitutes a huge technical challenge. Using the GPU has led to a very significant reduction of the registration time. The explanation of this reduction in processing time must be found in the parallelized architecture of the GPU.

An acceleration in the magnitude presented here is not only possible for the the Horn & Schunck method but should be attainable for other registration methods which lend themselves to being split into a large number of independent calculations. We expect the registration time of the methods presented in this paper to be reduced even more as the performance of graphics hardware increases. In the current implementation a fixed number of computations is performed for each image resolution. We plan to adjust the amount of computation to the degree of deformation actually present in the input images. This should also reduce computation time for most images. By splitting the images into blocks to be registered, is it also possible to distribute computations onto multiple GPUs.

This does however introduce an overhead from the memory synchronization needed at block boundaries so whether or not this would speed up the registration significantly is unclear.

So far we have shown that using GPUs the mentioned registration methods can be accelerated to a level which is acceptable for use in an online setting in which the deformable registration is done while the patient is still on the treatment couch. This is the first step towards online dose plan adjustment.

The processing power of GPUs can be utilized not only for registration as presented here, but for many of the compute intensive imaging tasks in IGRT making it an ideal and cost-efficient tool, which can help us getting further towards online IGRT.

Figure 1 :

 1 Figure 1: Visualization of the positioning of the landmark points on the images acquired from the head and neck cancer patient.

Figure 2 :

 2 Figure 2: The effect of applying the GREYCstoration filter to a head and neck CBCT image. The unfiltered image is the one on the right. The level of noise has been reduced without blurring the image. Window-level settings have been set to emphasize the difference between the two images.

Figure 3 :

 3 Figure 3: Differences between source and reference image before (left) and after (right) registration of an image from the POPI data set. A saggital slice and an axial slice are shown before and after registration. The source image is shown in a reddish color while the reference image is shown in a bluish color giving a gray scale image where the images align.

Figure 4 :

 4 Figure 4: Red/blue visualization of the difference between the rigid registration (left) and the deformable registration (right) of CBCT image 3 to CBCT image 1. A saggital slice and an axial slice are shown for each registration.

Figure 5 :

 5 Figure 5: Illustration of the transformation applied to the axial slice shown in figure 4.The transform is used to deform a rectilinear grid with a grid spacing of 10 mm. Only the in-plane deformation can be seen.

Figure 6 :

 6 Figure 6: Red/blue visualization of the difference between the rigid registration (left) and the deformable registration (right) of CBCT image 6 to the planning CT image. A saggital slice and an axial slice are shown for each registration.

Table 1 :

 1 Target registration error (TRE) compared to original distances of landmark points in the POPI data set.

			Mean Orignal
	Registration		Distance /		Mean TRE /
	(source/reference	Standard deviation	Maximum Orignal	Standard deviation	Maximum TRE
	image number)	(mm)	Distance (mm)	(mm)		(mm)
	0/1		0.5 / 0.5		2.4	0.7 / 0.3	1.4
	2/1		0.5 / 0.6		2.6	0.7 / 0.4	1.7
	3/1		2.2 / 1.8		6.6	1.3 / 0.8	3.4
	4/1		4.5 / 2.5		10.0	1.2 / 0.5	2.7
	5/1		6.0 / 2.9		12.1	1.3 / 0.7	3.6
	6/1		6.5 / 3.3		14.0	1.2 / 0.6	3.2
	7/1		5.5 / 3.0		14.0	1.3 / 0.6	2.8
	8/1		3.8 / 1.6		6.2	1.0 / 0.5	2.7
	9/1		2.1 / 1.0		4.5	0.9 / 0.6	2.5
	Average		3.5 / 2.0		8.0	1.1 / 0.6	2.7
							After deformable
			Before registration	After rigid registration	registration
	(source/refer	Mean / Std.		Mean / Std.		Mean / Std.
	ence image)	dev. (mm)	Max (mm)	dev. (mm)	Max (mm)	dev. (mm)	Max (mm)
	2/1	4.3 / 1.3	6.4	1.7 / 1.1	3.5	1.6 / 0.4	1.9
	3/1	4.2 / 1.6	6.5	2.5 / 1.1	3.7	1.8 / 1.0	3.4
	4/1	5.2 / 1.1	6.6	1.8 / 1.2	3.9	1.4 / 1.0	3.1
	5/1	6.5 / 0.8	8.0	1.4 / 0.9	2.7	1.4 / 0.6	2.1
	6/1	7.2 / 0.9	8.6	1.7 / 0.7	2.8	1.7 / 0.8	3.0
	Average	5.8 / 1.1	7.2	1.8 / 1.0	3.3	1.6 / 0.8	2.7

Table 2 :

 2 Evaluation of distances of bony landmarks for the CBCT-to-CBCT registration.

		After rigid registration	After deformable registration
		Mean / Std. dev.		Mean / Std. dev.
	Source image	(mm)	Max (mm)	(mm)	Max (mm)
	CBCT 1	2.2 / 1.0	3.7	1.4 / 0.6	2.7
	CBCT 2	1.7 / 0.4	2.2	1.4 / 0.5	1.9
	CBCT 3	1.9 / 0.5	2.9	1.7 / 0.3	2.2
	CBCT 4	2.8 / 0.6	4.1	2.7 / 1.1	4.9
	CBCT 5	2.0 / 0.3	2.5	1.7 / 0.5	2.6
	CBCT 6	2.3 / 0.7	3.7	2.0 / 0.6	3.3
	Average	2.2 / 0.6	3.2	1.8 / 0.6	3.0

Table 3 :

 3 Evaluation of distances of bony landmarks for the CBCT-to-CT registration.

Acknowledgements

The work presented in this paper has been partly funded by Varian Medical Systems, Inc., Palo Alto, California. We would like to thank the anonymous reviewers for helping us improve this paper.