Benoît Larrat
email: benoit.larrat@espci.fr

Mathieu Pernot

Gabriel Montaldo

Mathias Fink

Mickaël Tanter

MR-guided adaptive focusing of ultrasound

Keywords: Adaptive focusing, Hadamard, MR Acoustic Radiation Force Imaging, MR guided HIFU, phase aberrations

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Adaptive focusing of ultrasonic waves is a challenging problem in the field of medical ultrasound. The heterogeneities of biological tissues in terms of speed of sound, density or ultrasonic absorption induce a distortion of the ultrasonic wavefield that can result in a partial destruction of the focusing pattern. In order to restore the focusing quality, adaptive focusing relies on the use of ultrasonic arrays to correct the distortions induced by the propagation medium. This correction is performed by estimating and applying different time shifts (or phase shifts for monochromatic waves) on each element of the array. In some situations, when the ultrasonic array relies both on transmit and receive channels, it may be possible to rely on the echoes of a bright reflector or a point like active source located inside the biological tissues. From the signals received on the array, the so called Green"s function, one has to time reverse the wavefield (or phase conjugate for monochromatic signals) to focus back on the initial position. For example, time reversal focusing [START_REF] Fink | Time reversal of ultrasonic fields. I. Basic Principles[END_REF], [START_REF] Wu | Time reversal of ultrasonic fields. Il. Experimental results[END_REF] has encountered a significant success in the ultrasound community because piezoelectric transducers enable the exact recording and transmitting of broadband signals (i.e. transmission and recording of both amplitude and phase of the wavefield). Time reversal focusing was also extended to heterogeneous and dissipative media by applying the concept of spatio-temporal inverse filter [START_REF] Tanter | Time reversal and the inverse filter[END_REF]- [START_REF] Aubry | Optimal focusing by spatiotemporal inverse filter. II. Experiments. Application to focusing through absorbing and reverberating media[END_REF].

In a weak scattering medium with no strong reflector which is the case in many biological tissues, it is also possible to take advantage of the spatial coherence of the backscattered signals to estimate the time shifts induced by a thin aberrator with slow spatial variations [START_REF] Donnell | Phase aberration measurements in medical ultrasound: Human studies[END_REF]- [START_REF] Astheimer | Reduction of variance in spectral estimates for correction of ultrasonic aberration[END_REF].

However, the direct measurement of the Green"s function between the target location and the ultrasonic array is not always possible. In some configurations, echoes backscattered by the medium heterogeneities are too weak to be recorded by the array. For example, focusing ultrasound in the brain for therapeutic or imaging purposes has long been limited by the strong heterogeneities of the skull bone [START_REF] Fry | Acoustical properties of the human skull[END_REF]. The large discrepancy between the high speed of sound of the skull (up to 4000m/s) and the low speed of sound of brain tissues (about 1500 m/s) completely destroys the focusing pattern in the brain. The strong attenuation from the skull also makes backscattered signals from brain tissues unusable [START_REF] Pernot | High power transcranial beam steering for ultrasonic brain therapy[END_REF]. Similar issues are found, to a less extent, in other deep organs located under strong aberrating tissue layers such as the liver which is partially covered by the ribs [START_REF] Aubry | Transcostal high-intensityfocused ultrasound: adaptive focusing feasibility study[END_REF], [START_REF] Ballard | Adaptive Transthoracic Refocusing of Dual-Mode Ultrasound Arrays[END_REF]. In other configurations, the ultrasonic array can only be used in transmission and no receive channels are available to record backscattered echoes. This is very often the case in High Intensity Focused Ultrasound (HIFU) where the therapeutic array is only devoted to therapy and consequently is driven by electronic channels working only in the transmit mode. One can no more rely on ultrasonic echoes carrying information about the aberrations of the propagation medium.

For such configurations, we have recently introduced a novel method called energy-based adaptive focusing [START_REF] Herbert | Energy based adaptive focusing: application to non invasive aberration correction of ultrasonic wavefield[END_REF]. The general principle, which can be applied to any kind of waves in physics, relies on the indirect estimation of the wave intensity at the target for different coded excitations in order to obtain the time shifts (or phase shifts) information correcting the aberrations. By transmitting Hadamard-coded signals with an array of transducers and estimating the beam intensity at the target, this approach was shown to achieve a direct and accurate phase aberration correction without any phase measurement. In medical ultrasound, this idea is particularly interesting because ultrasonic waves interact with biological tissues through physical effects linked to the wave intensity such as the acoustic radiation force or tissue heating due to the absorption of ultrasound. Thus, the quantitative measurement of tissue displacement or temperature elevation [START_REF] De Senneville | Magnetic resonance temperature imaging[END_REF] at the target can be used for the indirect estimation of the local beam intensity.

The use of tissue displacements induced by the acoustic radiation force of the ultrasonic beams for energy-based adaptive focusing was validated by Herbert et al [START_REF] Herbert | Energy based adaptive focusing: application to non invasive aberration correction of ultrasonic wavefield[END_REF]. Due to the transfer of momentum from the acoustic wave to the medium [START_REF] Torr | The acoustic radiation force[END_REF], the radiation force induced by a short burst of ultrasound (~100µs) generates a local tissue displacement of several tens of micrometers at the focus of the ultrasound beam [START_REF] Sarvazyan | Shear wave elasticity imaging: a new ultrasonic technology of medical diagnostics[END_REF]. Small displacements (typically tens of microns) can be easily detected by ultrasonic means while remaining under the safety norms for medical ultrasound imaging. The use of this acoustic radiation force as a source of shear vibrations is currently used in clinical applications for the estimation of local stiffness [START_REF] Nightingale | Acoustic radiation force impulse imaging: In vivo demonstration of clinical feasibility[END_REF]- [START_REF] Muller | Quantitative Viscoelasticity Mapping Of Human Liver Using Supersonic Shear Imaging: Preliminary In Vivo Feasability Study[END_REF]. In our previous study [START_REF] Herbert | Energy based adaptive focusing: application to non invasive aberration correction of ultrasonic wavefield[END_REF], tissue displacement was estimated using speckle tracking obtained from pulse echo ultrasound sequences which provided quantitative and real time measurements of tissue displacement. The ultrasonic array was able both to transmit the ultrasonic coded excitations and to estimate the resulting displacements at the targeted location for a complete auto-correction of the aberrations. However, although ultrasound-based motion estimation was possible in this "proof of concept" experiment, ultrasound imaging cannot be used accurately through very strong phase aberrating medium. Indeed, beam distortions dramatically degrade the efficacy of speckle tracking at the target location. Moreover, in complex cases such as transcranial focusing, one cannot even rely on ultrasonic backscattered echoes due to the strong attenuation of the skull bone. In such cases, another imaging modality is required to provide an indirect estimation of the beam energy at the target.

We propose here to use Magnetic Resonance (MR) imaging to guide the energy-based adaptive focusing of the ultrasonic beam. Although conventional MR systems cannot image the very small displacement (on the order of 1 nm) induced by ultrasound at high frequency (~ 1MHz), MR imaging can be used to map larger displacement (on the order of 1µm) induced by the acoustic radiation force of ultrasonic longer bursts (typically some hundreds of microseconds) at a much lower frequency (< 1kHz). In recent studies, MR motion sensitive sequences have been developed to map displacements induced by the acoustic radiation force in biological tissues [START_REF] Sarvazyan | Shear wave elasticity imaging: a new ultrasonic technology of medical diagnostics[END_REF], [START_REF] Sinkus | Potential of MRI and Ultrasound Radiation Force in Elastography: Applications to Diagnosis and Therapy[END_REF]- [START_REF] Mcdannold | Magnetic resonance acoustic radiation force imaging[END_REF]. By synchronizing the short ultrasound bursts at precise timings of the MR sequence it is possible to detect small tissue displacements on the order of one micrometer with millimeter spatial resolution. To our knowledge, the experiments presented in this paper correspond to the first demonstration of MRI adaptive focusing of ultrasonic beams. This approach is particularly interesting for HIFU applications with therapeutic arrays used only in transmit mode as many HIFU therapeutic applications are already being performed under MR guidance, such as uterine fibroids [START_REF] Fennessy | A review of magnetic resonance imaging-guided focused ultrasound surgery of uterine fibroids[END_REF], breast tumor [START_REF] Gombos | Breast focused ultrasound surgery with magnetic resonance guidance[END_REF] or brain tumor ablation [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]. Although not widely available clinically yet, these treatments currently use MR to target the ablation zone and monitor the temperature elevation. In this paper, we demonstrate that MR guidance of HIFU therapy can also be used to achieve the adaptive focusing of the ultrasound beam prior to the treatment. This technique could potentially improve the focusing quality in all HIFU applications [START_REF] Kennedy | High intensity focused ultrasound: surgery of the future?[END_REF]. In the case of transcranial brain therapy, this technique could have even more importance. In practical clinical configurations, it could potentially become a very efficient and fast way to provide non invasive adaptive focusing of the ultrasonic beam through the skull, just before the treatment.

A theoretical framework of the energy-based adaptive focusing is first developed in the Material and Methods section and the link between energy based and time reversal focusing is illustrated. The experimental setup and the MR motion sensitive sequences are then presented. In the second section, the results of the phase aberration corrections obtained through strong aberrating layers are shown. The performance of the aberration correction process is investigated as a function of the motion encoding duration and of the number of available independent channels. Finally, the limitations and potential applications are discussed in the last section.

Materials and Methods

A -Theory

As explained above, MR systems are not sensitive enough and do not provide a sufficient temporal resolution to detect the very small displacements (nm) at ultrasonic frequencies (MHz) but can detect larger displacements (some µm) as those generated by the radiation pressure. By focusing a long burst of ultrasound focused at a given location, it is possible to create a volumic radiation force inside a dissipative medium at relatively low frequency (low kHz range). This force is due to the momentum transfer from the acoustic wave to the medium, caused by dissipation or reflection of the acoustic wave [START_REF] Torr | The acoustic radiation force[END_REF]. In medical applications, dissipation remains the dominant term in the radiation force and micrometric displacements can easily be induced in tissues with limited acoustic intensities.

For a dissipative medium, it is common to write that the radiation pressure is a volumetric force proportional to the acoustic intensity:

e I c e c P F       2 2 2 0 2   , (1)
where P , c and are respectively the pressure, the velocity and the attenuation of ultrasound,  being the mass density of the medium. e  is the direction of ultrasonic propagation and

) 2 /( 2 c P I  
is the local acoustic intensity. This force will generate bulk and shear waves propagating in the medium with their respective speeds, depending on the elastic properties of the medium. The spatio-temporal distribution of the displacement field generated by this force is complex and depends on many parameters such as the size and shape of the radiating focal spot and the viscoelastic properties of the medium [START_REF] Bercoff | The role of viscosity in the impulse diffraction field of elastic waves induced by the acoustic radiation force[END_REF]. The maximum displacement is found to be linearly linked to the intensity.

As monochromatic signals are considered here, it is convenient for sake of clarity to use the Fourier transform of the signals. Indeed, in the Fourier domain, the wave propagation between the ultrasonic array and the focal plane can be described very simply in a matrix formalism. For more details, the reader can refer to [START_REF] Tanter | Time reversal and the inverse filter[END_REF]. Working with monochromatic waves, the transmit signal of each transducer X i (1 i N) can be expressed by a complex number i e describing the amplitude and the phase of the emission. The total pressure received at the desired focus can be expressed as a complex number p containing the phase and amplitude of the local monochromatic field:

   N i i i e g p ... 1 , (2)
where i g are complex coefficients that relate the emission signal e i on each element X i with its resulting field at the focal point. In the Fourier domain, these coefficients are formally describing the Green"s function between the positions of the transducers i x and focal point , F r g i =Gr F ,x i ().

In a compact notation, an emission vector (

) 3
In the proposed technique, the pressure at the focus is maximized by adjusting the phases and amplitudes of the emission signals. This maximization has to be done using only the information on the intensity (i.e. the displacement) at the focal point.

The most straightforward method [START_REF] Urban | Phase aberration correction using ultrasound radiation force and vibrometry optimization[END_REF] consists in a two-transducer optimization. If one considers a configuration where the first transducer (working as a reference transducer) transmits simultaneously with another transducer that is phase shifted relatively to the reference, the emission

vector corresponds to ] 0 .. 0 , , 0 .. 0 , 1 [  j e  e
. The pressure fields generated by both transducers interfere at the desired focal point. It is possible to tune the phase  until the optimal phase is found which maximize the intensity in the desired focal point. The procedure must be repeated with each transducer of the array.

However, the radiation force generated by two single elements of the array is not sufficient to induce a displacement detectable by the MRI system. To solve this problem, Herbert et al. proposed to apply the same optimization principle to large sets of elements, the so-called "virtual transducers", instead of applying it to two single elements. The phase shift between the emissions of those two virtual transducers is then optimized. In practice, all elements can be used to build different virtual transducers. Indeed, an elegant set of independent emission vectors is given by the Hadamard"s matrix,

] ..h h , [h H n 1 2



, where each column i h of the matrix represents an emission vector (see Figure 1). The pressure obtained at the focus from the emission of

g.H p  H . ( 4
)
The goal of the aberration process is to estimate the vector g that describes the propagation from the array to the focal point. The knowledge of g will provide us a direct way to correct for the aberrations. For this purpose, the phase of each emission set i h has to be optimized relatively to the reference virtual transducer 1 h . To do that, the intensity resulting of the interferences between the reference emission 1 h and each of the others emission signals i h is measured. A practical method to recover the amplitude and phase of

H i
p is to perform the following four particular combined emissions:

) ( ) ( ) ( ) ( 1 2 2 1 2 2 1 2 1 1 2 1 i d i i c i i b i i a i j j h h s h h s h h s h h s         . ( 5
)
These emissions present the advantage of giving enough displacement to be measured by the MRI system because almost half of the transducers of the array are used simultaneously. Another interest is that the relative amplitude level required from each transducer is always 0 or 1 (see example in fig. 1) and the electronic channel does only control the relative phase between the transducers. This simplification has an important impact in terms of electronics complexity. Indeed, this kind of transmit codes can be provided even by HIFU hardware systems that do not allow variable transmit voltage on each individual electronic channel. . So, intensities generated by the four different excitations described above are:

) /( ) ( ) /( ) ( ) /( ) ( ) /( ) ( ] [ ] [ ] [ ] [ c p p imag p p I c p p imag p p I c p p real p p I c p p real p p I H i H H i H d i H i H H i H c i H i H H i H b i H i H H i H a i                 2 2 2 2 1 2 2 1 2 1 2 1 1 2 2 1 2 1 2 1 1 2 1 2 4 1 2 1 4 1 1 2 1 2 4 1 2 1 4 1 . ( 6
)
The pressure field at the focal point is deduced from this system of equations:

) ( ) ( ) ( ) ( c i d i H i H b i a i H i H I I c p p imag I I c p p real       1 1 2 . ( 7
)
As H p 1 is the reference signal obtained at the focus, we can consider arbitrarily that the phase of H p 1 is 0 and its relative amplitude can be obtained as (2cI 1 a ) 1/2 by using

1 h h i 
in equation [START_REF] Donnell | Phase aberration measurements in medical ultrasound: Human studies[END_REF].

Finally, the different pressures received at focus are:

) /( ) ( ) ( / ) ( ) ( H c i d i H i H b i a i H i p I I p imag p I I p real 1 1 2     (8)
These values of H i p can be estimated and used in equation [START_REF] Tanter | Optimal focusing by spatiotemporal inverse filter. I. Basic principles[END_REF]. By inversing the Hadamard matrix (H - 1 =H t ), the propagation vector (both in amplitude and phase) describing the Green"s function between the array and the desired focal spot is finally:

t H .H p g  . ( 9
)
Once the Green"s function is estimated, it is very simple to correct for the aberrations. It was shown that the optimal way to maximize the ratio between the acoustic intensity received at focus and the total acoustic intensity delivered by the array is to use the time reversed version of the Green"s function [START_REF] Tanter | Time reversal and the inverse filter[END_REF]. This consists in inverting the phase and using the amplitude information of the Green"s function as a transmission signal. Indeed, time reversal focusing performs both a spatial and temporal matched filter for the propagation operator [START_REF] Tanter | Time reversal and the inverse filter[END_REF].

Regarding temporal focusing, this demonstration is straightforward and is recalled here. For a set of transmitted signals s, the ratio R between the acoustic intensity received at focus and the total acoustic intensity delivered by the array is given by :

2 2 2 2 2 g s gs s    p R . ( 10
)
Thanks to Cauchy-Schwartz relation, ratio R becomes maximal when the transmission signal s is proportional to the conjugate of the Green"s function:

i i F i G g X r s   ) , ( (11)
The same kind of arguments can be proposed both for space and time and the interested reader can refer to [START_REF] Tanter | Time reversal and the inverse filter[END_REF] for a complete demonstration. Applying an amplitude i g and a phase value -  i g  on each element i will maximize the ratio R and thus the efficiency of the acoustic deposit at the focus. It corresponds to the solution provided by the well known time reversal approach [START_REF] Tanter | Time reversal and the inverse filter[END_REF]. One can also decide to reduce the aberration correction to a single phase compensation (by applying the same amplitude on each element of the array and only inverting the phase of the Green"s function). This approach increases the amount of energy deposit at focus but could result in undesired heat deposit along the travel path in highly absorbing areas.

B -Experimental setup

Figure 2 shows the experimental setup. A 128-element linear phased array (Vermon, France) operating at 6MHz (spatial pitch 0.3 mm) was used inside a 7T MRI scanner (Bruker, Germany). The array is controlled by a fully programmable electronic system (Lecoeur Electronics, France) able to emit independent signals on 64 US channels. We thus used only one half of the 128-element array.

The ultrasonic emissions were directly triggered by the MR scanner itself. It allowed an accurate synchronization of the US emissions with the MR acquisition. The ultrasonic waves propagated inside a tissue mimicking phantom made of 4% gelatin. An aberrator was placed between the array and the phantom. Depending on the experiment, the phase aberration was either "virtual" using a programmable aberrator or "real" using a physical aberrator. The programmable aberrator consisted of additional programmable delays placed at the emission stage in the US electronics. The physical aberrator consisted of two stripes of PolyDiMethylSiloxane (PDMS) or various thickness and width placed between the array and the gelatin phantom. This material was chosen for its low attenuation and its speed of sound (approximately 1080 m/s) which is significantly lower than in the phantom. A burst of ultrasound of typically 400µs generates displacements in the medium through the radiation force effect. These displacements are relatively moderate (a few µm) and as soon as ultrasounds are switched off, the medium gets back to its original position while a radiated shear wave propagates in the elastic gel within a few ms.

C -Motion sensitive MR sequences

As recently published [START_REF] Sarvazyan | Shear wave elasticity imaging: a new ultrasonic technology of medical diagnostics[END_REF], [START_REF] Sinkus | Potential of MRI and Ultrasound Radiation Force in Elastography: Applications to Diagnosis and Therapy[END_REF]- [START_REF] Mcdannold | Magnetic resonance acoustic radiation force imaging[END_REF], ultrasound induced tissue displacements can be mapped via MRI. To image these displacements, the combined MRI/ultrasonics sequence must take into account two important points: it must be sensitive enough to detect a micrometric displacement and it must be fast enough to image it before shear wave propagation occurs. Shear wave propagation degrades the spatial localization of the detailed focusing pattern and makes mechanical energy appear where no radiation force is present. The sequence is schematized in Figure 3. A spin echo MR sequence was modified to encode for motion (TE=15ms/TR=500ms). This consisted in adding 2 periods of sinusoidal gradients symmetrically positioned before and after the 180° radiofrequency pulse of the sequence. A frequency of 600Hz was chosen for the gradients with maximum amplitude of 300mT/m. As detailed latter, this choice maximized the sensitivity to motion while minimizing the US emission duration. The displacements were mapped in one slice oriented transversely to the US beam. The slice was positioned to contain the chosen point where the aberration correction is to be performed. The motion sensitizing gradients (MSG) were placed on the slice selection axis to be sensitive to the displacement along the US beam. The radiation force induced displacement at the focus was assumed to be mostly oriented along that direction. The displacement field was mapped in a 32*256 voxel matrix with an in-plane resolution of 300µm*200µm and a slice thickness of 3mm. The sensitivity to motion was about one micron. Two 400µs ultrasound bursts were emitted at the beginning of the second half period of each MSG. Such a configuration avoids any motion during the first lobes of the MSG which would destroy motion encoding. We found it to be the best configuration both in term of motion encoding and in term of signal to noise ratio. The acoustic duty cycle remained below 0.2% which ensured the safety of the US transducer as well as the absence of heat deposition in the phantom. The MR phase ΔΦ is linked to the spatiotemporal displacement profile u(r,t) via the following relationship:

, ( 12
)
where  is the proton gyromagnetic constant and G the motion sensitizing gradient magnetic field.

D -Adaptive focusing process

In this work, only 64 adjacent channels were used to keep the acquisition time reasonable.

Furthermore, the focal plane was located around 2.5cm away from the array. Given the directivity of the transducers and the pitch of the array, very little signal would be emitted from the lateral elements if more than 64 channels were used and thus the phase correction would become difficult for those elements.

The three delays described in Figure 2 enabled to decompose the total delay applied to each channel during the acquisitions. Starting with a plane wave, a prefocusing cylindrical delay was added to target one chosen point in space where the correction is to be performed. This prefocusing delay generated a sharp focus at the target. Then, as shown in the theoretical part the adaptive focusing procedure implied to emit coded emissions. This consisted of using only some of the channels with additional phase delays for some of them according to equation 5. Then a phase aberrator was either programmed (see section Results -C) and added numerically to the previous delays or physically placed between the ultrasonic array and the gelatine medium (see section Results -D).

The excitation-imaging sequence (see Figure 3) was repeated 4*64 times corresponding to the four emissions (s a , s b , s c and s d ) in equation ( 5) for the 64 Hadamard vectors. The emission amplitude was set to the maximum value deliverable by the electronics for all acquisitions (80V peak to peak).

For the 256 emissions, the whole acquisition lasted around 2.5 hours. The obtained MR phase maps were unwrapped. A reference image acquired without ultrasound was subtracted to each image to get real displacement weighted data and quantify the acoustic power at the focus. The set of time delays and amplitudes to be applied on each element of the US array for the optimal focusing was then recovered by the direct inversion of equation 9. It is to be noticed that the MR phase was recorded over a 2D plane while the sole knowledge of the acoustic intensity at the targeted point was really needed. The 2D imaging capability, although time consuming could be used in different manners as discussed later.

In a first experiment, an aberrating phase profile with strong phase shifts was programmed in the transmitted US beam. The recovered phase correction was compared to the programmed input. In a second experiment, a physical aberrator was placed between the array and the phantom, in this case, the recovered phase correction was compared to the Green"s function of the targeted point as measured via a pulse echo experiment with an acoustic reflector placed at the target.

E -Simulations

Simulations of the elastic response of the gel were also performed in order to determine the influence of the duration of motion encoding on the phase correction. The theoretical visco-elastic

Green"s function was computed using the method developed by Bercoff et al. [START_REF] Bercoff | The role of viscosity in the impulse diffraction field of elastic waves induced by the acoustic radiation force[END_REF]. The freely available ultrasound simulation program Field II [START_REF] Jensen | Calculation of pressure fields from arbitrarily shaped, apodized, and excited ultrasound transducers[END_REF]- [START_REF] Jensen | Field: A Program for Simulating Ultrasound Systems[END_REF] was used to simulate the pressure field generated by the linear array in monochromatic regime (6MHz). The resulting acoustic radiation force was computed over a 4x4x4 mm3 grid with a spatial step of 0.25 mm.

Then, a dedicated code [START_REF] Bercoff | The role of viscosity in the impulse diffraction field of elastic waves induced by the acoustic radiation force[END_REF] was used to simulate the viscoelastic displacement response u(r,t)

of the gel submitted to a 400µs sonication. This method is based on an explicit derivation of the 3D

Green's function in a viscoelastic medium taking into account shear, bulk, and coupling waves. The acoustic radiation force was used as the volumic shear source. The shear wave speed of the medium was taken as 1.5m/s as measured on 4% gelatin gels via the Supersonic Shear Imaging technique [START_REF] Bercoff | Supersonic shear imaging: a new technique for soft tissue elasticity mapping[END_REF].

The shear viscosity was neglected. The mass density was taken as 1kg/m3 and the speed of sound 1540m/s.

Finally, the complete spatio-temporal viscoelastic displacement response u(r,t) was used to simulate the MR motion encoding. The MR phase was computed using equation 12 with sinusoidal MSG gradients of varying frequency. It allowed studying the influence of the speed of motion encoding (i.e. the frequency of the MSG) on the resulting recovered phase. The MR phase was then discretized over a grid composed of voxel of 0.25x0.25x3 mm3. The influence of the number of channels used for the acquisition was studied numerically.

Results

A -Experimental verification of the relation between MR signal and acoustic intensity

An important assumption of the energy based adaptive focusing method is that the measured signal at the targeted point is linearly linked to the acoustic intensity. This is verified through the recording of the viscoelastic response of the phantom for different amplitudes of the emitted ultrasonic signal. This relationship was measured prior to all experiments and the displacement was found to rise with the applied voltage U according to a quadratic law as shown on Figure 4 for an experiment with a physical aberrator in front of the array (see Results-C). From the integration over time in equation ( 12), the MR phase at the targeted point is linearly linked to the acoustic intensity I so that ΔΦ~U²~I.

On Figure 4, the relatively poor fit of the data at low emitted power is explained by the limited signal to noise ratio of the MR phase images. All subsequent acquisitions were performed at maximum power.

B -Simulations

In order to demonstrate the feasibility of the experiment and to determine the influence of the frequency of the MSG on the results, the whole experiment was simulated. The prefocusing delay was chosen to target a point located 30mm deep in the medium. The pressure field was simulated in 3D for each of the 256 emissions with 0.25mm*0.25mm*0.25mm resolution using the Field II software. The corresponding viscoelastic responses were obtained for a 400µs ultrasonic burst. The MR motion encoding step was simulated for 3 different sinusoidal MSG frequencies assuming that the US bursts were emitted at the beginning of the second half periods on Figure 3. Figure 5 shows the simulated MR phase maps obtained for two different ultrasonic emission phase profiles. As both the emission delays are Hadamard-coded transmissions and the ultrasonic beam propagates through the aberrator, the ultrasonic beam is not sharply focused at the target. The spatial distribution of the pressure field (and the acoustic radiation force) is complex and presents several important lobes as shown on the right column of Figure 5. Figure 5 shows that for some emission phase profiles such as emission #36, the spatial peak of the radiation force distribution is not even located at the target. In order to estimate the correct beam intensity at the target, the tissue displacement should reproduce well enough the spatial distribution of the acoustic radiation force. This assumption is correct only during a very short time at the beginning of the mechanical excitation. However, for long excitations, the spatial distribution of tissue displacement evolves rapidly due to the viscoelastic response of the medium. give very good results and was adopted for the following experiments. In theory an even shorter duration should further improve the results for a 400µs burst assuming an infinite signal to noise ratio on the MRI acquisitions. However, in 4% gelatin gels and given the available acoustic power, 600Hz appeared to be the experimental limit to get a reliably quantifiable MR phase image. In general, the choice of the MSG frequency remains a compromise between shear wave propagation and MR motion sensitivity.

Finally, the focusing pattern recovered after adaptive focusing correction for different encoding speed was compared to the ideal focusing pattern. As shown on figure 7, the simulation at 600Hz is very close to the ideal case whereas the focus quality is strongly destroyed with increasing encoding duration. are compared with ideally corrected (green curve) and not corrected (blue curve) ones.

From those simulations, a frequency of 600Hz was chosen as a good compromise between encoding speed and signal to noise ratio.

C -Experiments with a programmable aberrator

Using programmable delays (see Figure 2) the presence of a known phase aberrator between the array and the propagating medium is simulated. The difference with a real physical aberrator is that the phase aberration layer is infinitely thin and has no attenuation. MR imaging allowed accurate estimation of the acoustic power at the focus of the array. Two examples of MR phase maps recorded for the Hadamard coded emissions number 36 and 99 are given on Figure 8. After the complete 256 transmit experiment, the proposed adaptive focusing method was able to recover the spatial distribution of phase aberrations in gelatin phantoms for the same numerical phase aberrating layer as the one previously used for simulations. The programmable phase aberration profile is shown together with the recovered phases on Figure 9. A low error was found on the non invasive phase aberration estimation (standard deviation of 0.30 radians). A refocusing test was done after phase correction to compare the technique with an ideally focused beam. An estimate of the acoustic power at the focus was mapped using MR with and without phase correction and an increase of a factor of 10 was found. Figure 10 shows the measured MR phase along the X axis for the non aberrated focusing phase profile, the aberrated phase profile and the corrected phase profile. It is to be noticed that the focal spot after phase correction had a slightly greater energy than with a perfect non aberrated focusing law. This may be caused by the unexpected presence of a small physical additional phase aberrations in front of the array (bubbles, impurities, transducer defects) which were also corrected by the adaptive focusing procedure.

D -Experiments with a physical Aberrator

The same procedure was employed with two PDMS stripes placed against the linear array (see Figure 11). The aberration profile recovered along the 64 elements after phase correction is given on Figure 12. The effect of the two stripes of different widths is clearly visible. To evaluate the actual phase aberration, the whole setup was taken outside of the magnet and a strong reflector (thin iron bar) was inserted in the gel at the approximate targeted point. A pulse echo ultrasonic experiment gave the Green"s function of the target on the 64 elements used for the MR guided adaptive focusing experiment. This aberration is also plotted on figure 12. A good accordance was found between both measurements (standard deviation of the difference equal to 0.83 radians). The procedure also gives a calculated profile for the amplitude aberration. The calculated correction for the amplitude shows very weak energy transmission in the region of the strongest aberration (elements on the right side of the array). This amplitude aberration explains the relatively high error on the estimation of the phase correction in that area. The proposed adaptive focusing procedure is capable of recovering the optimal amplitude to emit on each transducer in order to avoid absorbing and aberrating structures and maximize the energy at the focus. In this study, however, no amplitude correction was performed and the maximum emitted power was maintained identical on all 64 channels. Only the phase was corrected.

In order to assess the quality of the correction, the MR phase signal at the focus was measured after phase correction and an increase of a factor of 5 was found. Side lobes have disappeared. Figure 13 shows the measured MR phase in the focal plane along the direction of the array (X) for the aberrated phase profile and the corrected phase profile. Figure 14 shows the corresponding MR phase maps in the focal plane (XY).

E -Influence of the number of channels

The performance of the adaptive focusing method was investigated as a function of the number of independent channels. The complete adaptive focusing process was simulated with the previous programmable phase aberrator. The whole dataset corresponding to the 64 transducers elements (i.e. 4*64 simulated MR phase images) was first considered. Then, this dataset was truncated to select the 4*32 emissions corresponding to an adaptive procedure performed with only 32 transducers elements of twice the width of the previous individual elements. They corresponded to the Hadamard vectors for which 2 neighboring elements of the array always emitted the same signals.

Finally, in the same way, a second truncation was performed to simulate the procedure for N=16 transducer elements. In this case, the 64 elements were coupled by groups of 4 consecutive elements.

In this way, the total aperture of the array did not change but the spatial sampling was degraded. This means that for this particular aberrator, the optimal result could be reached with only 4*32 independent coded excitations as they will be sufficient to describe the full aberrations. This is of particular importance as the adaptive focusing procedure would require half the full acquisition time. The reduction from 64 elements to 32 elements of doubled size corresponds to a reduction of the inversion process to a lower rank in terms of independent spatial frequencies. As the results were unchanged here, the particular aberrator chosen for the experiment does not contain very high spatial frequencies. In other words, the correlation length characterizing the spatial distribution of the aberrations is higher than the doubled element size. This result clearly emphasizes the interest of studying the spatial frequency content of the aberrating medium along the array aperture. The previous example implies that the knowledge of the spectrum of a given aberration on the Hadamard"s vectors basis determines the number of independent spatial frequencies needed for an optimal correction. In this context, the use of the Walsh decomposition of the aberration is highly informative in order to optimize the experimental acquisition. Indeed, Walsh"s vectors are equal to

Hadamard"s vectors but are ordered from low frequencies to highly fluctuating signals. The decomposition of a given aberration profile on the Walsh basis gives its spatial frequency content as the first vectors of the Walsh basis correspond to low spatial frequency profiles and the last ones correspond to high frequency spatial fluctuations. This information is valuable because it enables to determine the optimal number of degree of freedom required in the array aperture to fully correct the aberrations.

The Walsh "spectrum" of the proposed numerical aberrator was calculated for the 64 channels and the degraded 32 and 16 channels simulated experiments (see Figure 16). As expected, it showed no peak over Walsh"s vector number 32 for the input aberration. This gives the minimum number of independent channels needed for the adaptive focusing in this special case. To evaluate the effects of the degraded acquisition, the recovered corrected focal spots using the phase corrections obtained with 64, 32 and 16 channels were compared (Figure 15.b). On the same plot, the simulated MR phase profile in the focal plane is given for the cases without aberrations and without any correction. The curves for 64 and 32 channels are strictly identical.

It is interesting to notice that the recovered phase profile for the 16 channels case has a Walsh spectrum which is more rich than a simple low pass truncation of the spectrum of the complete aberration (Fcut=16). The proposed physically based adaptive focusing procedure retrieves the best possible phase correction given the limited available information. This is suggested when comparing the energy at the focus for the 16 channels decimated acquisition with the one obtained from a theoretical phase correction given by a low pass truncation of the spectrum of the aberrator. The decimated acquisition gives better results.

Discussion

In this work, a new adaptive focusing technique is employed to correct the phase distortions of an ultrasonic wave beam. It is based on the measurement of the acoustic radiation force via MRI. The technique is successfully employed in simulations and real experiments for both programmable and physical aberrations with phase shifts up to 2π. The restored energy at the targeted point is measured and the maximum acoustic intensity is reached after correction. A sharp focus is restored and the intensity is increased by a factor of 5. Although the proof of concept of MR-guided adaptive ultrasonic focusing is illustrated using a 1D linear array performing the aberration correction in the 2D imaging plane of the ultrasonic array, this concept is straightforwardly applicable to any kind of 2D ultrasonic arrays.

The simulation work demonstrates the importance of motion encoding on the final result: due to shear wave propagation, the US bursts as well as the MR motion encoding must be as short as possible to keep an accurate mapping of the radiation force at the target. The stiffer the medium, the more significant this effect becomes. The experimental limitation comes from the MR sensitivity to motion which is related to the signal to noise ratio and to the motion encoding gradients of the sequence. At 7T, in relatively soft gelatin gels, the compromise was found for 400µs bursts and 600Hz

MSG. In the perspective of applying the proposed technique to transcranial adaptive focusing for brain HIFU, a gel with an elasticity of about 1.5m/s was employed. This value roughly corresponds to the published value in the brain over the same range of excitation frequencies [START_REF] Nicolle | Shear linear behavior of brain tissue over a large frequency range[END_REF].

In the present work, we reached a very good correction with a limited US duty cycle. Short bursts of 400µs were transmitted during the MR imaging sequence and the US duty cycle was kept below 0.2 % during the full acquisition time of 150 minutes. With such a low duty cycle, the temperature elevation at the end of the procedure was not significant which ensures the safety of this method.

The precision of the phase correction depends on the SNR of the displacement detection for all the coded transmissions. The retrieved delays can become erroneous for those of the elements which are much attenuated (amplitude aberration). This limit must be evaluated for each particular application. It has to be noticed that in such cases, the technique also finds that low amplitude signals have to be emitted by those attenuated channels.

This new adaptive focusing method should be applicable to clinical HIFU treatments in a whole body MR magnet. A major application of the proposed MR-guided approach should be transcranial brain HIFU. In this perspective, Table 1 sums up the values of the main technical parameters in a clinical system as compared with the one used in the present work. The ultrasound frequency is an important parameter to consider. Because the lateral dimension of the focal zone and the distance between the side lobes is the same order of magnitude than the ultrasound wavelength, the visco-elastic response of the medium requires more time to destroy the initial displacement pattern.

Therefore, the motion encoding duration can be increased when the US frequency decreases without significant blurring of the displacement maps due to wave propagation. This roughly compensates the lower available MSG strength at 1.5T and enables to reach similar sensitivity to motion in both cases.

The use of lower gradient strength for motion encoding is required for in vivo application considering MR safety guidelines [START_REF]Medical Electrical Equipment -Part 2-33: Particular requirements for the safety of magnetic resonance equipment for medical diagnosis[END_REF]. Moreover, the pressure amplitude required to induce micrometric displacements through the human skull bone at typical HIFU treatment frequencies (around 1 MHz) is of the order of some MPa and corresponds to typical acoustic intensities used in transcranial HIFU [START_REF] Pernot | High power transcranial beam steering for ultrasonic brain therapy[END_REF], [START_REF] Pernot | In vivo transcranial brain surgery with an ultrasonic time reversal mirror[END_REF][START_REF] Tanter | Compensating for bone interfaces and respiratory motion in high-intensity focused ultrasound[END_REF]. Moreover, the absence of skull displacement in the case of a stereotactic approach will ensure that the aberration correction can be performed without being degraded by motion artifacts.

Effects of the blood flow on displacement data are expected to be small and partly cancelled out by a systematic difference between Hadamard coded MR phase signals and zero emission signal.

In the framework of HIFU brain therapy, the aberrations introduced by the skull bone are similar than the ones corrected in the present experiments. The performances reached in the present study should thus be sufficient to correct aberrations from the skull in an acceptable manner. A common criterion states that 0.8 radians is an acceptable residual phase delays after correction [START_REF] Goodman | Statistical Optics[END_REF].

Here, we are below this limit in most cases. The acquisition time was not optimized in this work and the complete adaptive focusing sequence was lasting around 2.5 hours for 64 channels. This was mainly due to technical limitations of the electronics which enforced us to load at low transfer rate one US emission on the electronic boards at a time. This roughly doubled the total acquisition duration. Another limitation derives from the obligation to image a 2D plane with at least 32 lines with the available motion sensitive sequence.

Using a line scan approach is feasible with an orthogonal orientation for the excited and the refocused slices in the spin echo scheme. This will result in a significant time reduction of a factor 32. Shorter TR may also be used, especially in the brain where the MR signal is good. Last, 3N emissions should be theoretically sufficient to perform the proposed technique and recover the phase aberrator without loss in precision. All things considered, a time reduction factor of roughly 300 can be reasonably envisaged resulting in a 4 minutes acquisition time for the full adaptive focusing with 512 channels in a clinical MRI. Such a short acquisition time is totally compatible with clinical therapeutic transcranial ultrasound applications. In particular, time constrained applications become possible. For instance, ultrasonic treatment of stroke that requires a fast therapeutic response (<4 hours after the event) could be envisioned with such an adaptive focusing procedure [START_REF] Molina | Transcranial ultrasound in clinical sonothrombolysis (TUCSON) trial[END_REF].

In addition, the 2D imaging capability, although time consuming at first sight could be exploited to learn the correcting delays to properly focus at several different locations in parallel. For instance, it should be possible to focus at 2 different points in the same imaged plane with 2 sub-arrays resulting in 2 restored focusing capabilities, one with each half of the array. Even further, if the MRI is sensitive enough to motion, the array could emit a non pre-focused wave. After the transmission of all coded emissions, the knowledge of the acoustic energy at any specific point in the imaged plane should allow to find the right focusing law for this point.

Moreover, the Walsh analysis of the aberration gives another mean of optimizing the acquisition. For instance, it is possible to start the acquisition with the slow spatial phase fluctuations and determine a cutoff frequency over which no further correction is needed. This approach could help optimizing the acquisition time. It could also be used to perform a two step optimization with a rough low frequency correction followed by a refined high frequency one.

When some information is known on the spatial distribution profile of the aberrating layer, the adaptive focusing procedure can be accelerated as the number of required coded excitations is decreased. The example of transribs HIFU is a good illustration of this point. In transribs ultrasonic focusing, most of the aberration correction consists in shutting down the array elements that are located in front of each rib so that the ultrasonic beam only propagates in the intercostal spacing. As this intercostal spacing can be estimated before the treatment, the aberration profile could be potentially described by a limited number of Hadamard vectors and thus results in a short number of data acquisitions.

Finally, an important aspect of this adaptive focusing technique is related to safety issues in HIFU treatments. As the optimization process is based on the direct measurement of local displacements at the desired focus, it also provides a direct control of the efficacy of the adaptive beam correction process. The displacement field measured using the motion sensitive sequence can be used to check before the treatment that the correction is correctly maximizing the acoustic intensity delivered at the targeted focus.

1 

 1 describing the field received at focus for all the successive N emissions i h , one can write:

Figure 1 .

 1 Figure 1. Example of the Hadamard matrix for N=4 elements and the signals s i used for the optimization process with a total of 16 ultrasonic transmissions. The amplitude of s i is always 0 or 1.

Figure 2 .

 2 Figure 2. Experimental setup. Fully programmable electronics are able to send independant signals on

Figure 3 .

 3 Figure 3. Proposed MR sequence for radiation force imaging at 7T. Two 600Hz sinusoidal gradients

Figure 4 .

 4 Figure 4. Calibration curve obtained with a PDMS physical aberrator in front of the linear array. The

Figure 5

 5 Figure 5 compares three different MSG frequencies (with identical sonication durations). Only the

Figure 5 .

 5 Figure 5. Normalized simulated MR phase maps in XY plane with encoded displacement in the focal

Figure 6 .

 6 Figure 6. Numerical phase aberrations along the 64 elements linear array: input values (blue curve),

Figure 7 .

 7 Figure 7. Simulated profile of the focal spot in the X direction (along the array) measured with the

Figure 8 .

 8 Figure 8. Example MR phase images for two Hadamard coded emissions (#36 and #99) in the

Figure 9 .

 9 Figure 9. Numerical phase aberrations along the 64 elements linear array: programmable phase

Figure 10 .

 10 Figure 10. Retrieved MR displacement weighted signal at the target along the axis of the US array for

Figure 11 .

 11 Figure 11. Picture of the ultrasonic linear array with two PolyDiMethylSiloxane stripes positioned in

Figure 12 .

 12 Figure 12. Phase correction obtained with the proposed method for 2 physical aberrating stripes placed

Figure 13 .

 13 Figure 13. Measured MR phase in the focal plane along the direction of the array (X) for the non

Figure 14 .

 14 Figure 14. Measured MR phase maps in the focal plane (XY) prior any correction (left) and after the

Figure 15 .

 15 Figure 15.a presents the phase aberration profile recovered from the full 64 elements simulated

Figure 15 a

 15 Figure 15 a) Recovered phase aberration profiles from simulated experiments: input, recovered with

Figure 16 .

 16 Figure 16. Walsh spectrum of the programmed aberrator together with the spectra of the aberrators

Table 1 .

 1 Comparison of the main technical parameters in the perspective of a transposition of the technique at 1.5T.

	Magnetic field	7T	1.5T
	US probe	64 channels diagnostic linear	512 channels HIFU spherical
		array	array
	Medium	Gelatine 4%	Brain tissue
	US frequency	6MHz	1MHz
	MSG Gradients strength	300mT/m	40mT/m
	Encoding frequency	~600Hz	~100Hz