P Berceau

M Fouché

R Battesti

C Rizzo

Magnetic Linear Birefringence Measurements Using Pulsed Fields

In this paper we present the realization of further steps towards the measurement of the magnetic birefringence of the vacuum using pulsed fields. After describing our experiment, we report the calibration of our apparatus using nitrogen gas and we discuss the precision of our measurement giving a detailed error budget. Our best present vacuum upper limit is ∆n ≤ 5.0 × 10 -20 T -2 per 4 ms acquisition time. We finally discuss the improvements necessary to reach our final goal.

I. INTRODUCTION

Experiments on the propagation of light in a transverse magnetic field date from the beginning of the 20 th century. In 1901 Kerr [1] and in 1902 Majorana [2] discovered that linearly polarized light, propagating in a medium in the presence of a transverse magnetic field, acquires an ellipticity. In the following years, this linear magnetic birefringence was studied in detail by Cotton and Mouton [3] and it is known nowadays as the Cotton-Mouton effect. It corresponds to an index of refraction n for light polarized parallel to the magnetic field B that is different from the index of refraction n ⊥ for light polarized perpendicular to the magnetic field. For symmetry reasons, the difference between n and n ⊥ is proportional to B 2 . Thus, an incident linearly polarized light exits from the magnetic field region elliptically polarized. For a uniform B over an optical path L, the ellipticity is given by:

Ψ = π L λ ∆nB 2 sin 2θ, ( 1
)
where λ is the wavelength of light in vacuum, ∆n = nn ⊥ at B = 1 T and θ is the angle between light polarization and magnetic field. The Cotton-Mouton effect exists in any medium and quantum electrodynamics predicts that magnetic linear birefringence exists also in vacuum. It has been shown around 1970 [7,8] thanks to the effective Lagrangian established in 1935 and 1936 by Kochel, Euler and Heisenberg [4,5]. At the lowers two orders in α, the fine structure constant, ∆n can be written as:

∆n = 2 15
α 2 h3 m 4 e c 5 µ 0 1 + 25 4π α [T -2 ], ( 2
)
where h is the Planck constant over 2π, m e is the electron mass, c is the speed of light in vacuum, and µ 0 is the * Electronic address: carlo.rizzo@lncmi.cnrs.fr magnetic constant. The α 2 term is given in Ref. [7]. The α 3 term has been first reported in Ref. [9] and it corresponds to the lowest order radiative correction. Its value is about 1.5% of the α 2 term. Using the 2010 CODATA recommended values for the fundamental constants [14], Eq. ( 2) gives ∆n = (4.031699 ± 0.000002) × 10 -24 [T -2 ].

As we see, the error due to the uncertainty of fundamental constants is negligible compared to the error coming from the fact that only first order QED radiative correction has been calculated. The QED α 4 radiative correction should affect the fourth digit and the QED α 5 radiative correction the sixth digit. Thus, a measurement of ∆n up to a precision of a few ppm remains a pure QED test.

Experimentally, the measurement of the Cotton-Mouton effect is usually very challenging especially in dilute matter, thus all the more so in vacuum. Several groups have attempted to observe vacuum magnetic birefringence [12,13], but this very fundamental prediction still has not been experimentally confirmed.

Gas measurements date back to 1938 [6] and the first systematic work of Buckingham et al. was published in 1967 [10]. Investigations concerned benzene, hydrogen, nitrogen, nitrogen monoxide and oxygen at high pressures, and ethane. Since 1967, many more papers concerning the effect in gases have been published and Cotton-Mouton effect experiments have been employed as sensitive probes of the electromagnetic properties of molecules [6].

The measurement of the Cotton-Mouton effect in gases is not only important to test quantum chemical predictions. It is a crucial test for any apparatus which is dedicated to the search for vacuum magnetic birefringence. Measurement of the Cotton-Mouton effect in a gas is a milestone in the improvement of the sensitivity of such an apparatus. Typically measurements of the linear magnetic birefringence in nitrogen gas are used to calibrate a setup [12,13,15].

In the following we present magnetic linear birefringence measurements performed in the framework of our "Biréfringence Magnétique du Vide" (BMV) project. It is based on the use of strong pulsed magnetic fields, which is a novelty as far as linear magnetic birefringence is concerned, and on a very high finesse Fabry-Perot cavity to increase the effect to be measured by trapping the light in the magnetic field region. The use of pulsed fields for such a kind of measurements has been first proposed in Ref. [11]. In principle, pulsed magnetic fields can be as high as several tens of Tesla, which increases the signal, and they are rapidly modulated which decreases the 1/fflicker noise resulting in an increase of the signal to noise ratio. Both advantages are supposed to compensate the loss of duty cycle since only few pulses per hour are possible. A feasibility study, which discusses most of the technical issues related to the use of pulsed fields coupled to precision optics for magnetic linear birefringence measurements, can be found in Ref. [17].

In this paper we present the realization of further steps towards the measurement of the magnetic birefringence of the vacuum using pulsed fields. After describing our BMV experiment, we report the calibration of our apparatus with nitrogen gas and we discuss the precision of our measurement giving a detailed error budget. Finally, present vacuum upper limit is reported and we discuss the perspectives to reach our final goal.

II. EXPERIMENTAL SETUP AND SIGNAL ANALYSIS

A. Apparatus

The BMV experiment is detailed in Ref. [17]. Briefly, as shown on Fig. 1, 30 mW of a linearly polarized Nd:YAG laser beam (λ = 1064 nm) is injected into a Fabry-Perot cavity consisting of the mirrors M 1 and M 2 . The laser frequency is locked to the cavity resonance frequency using the Pound-Drever-Hall method [18]. To this end, the laser is phase-modulated at 10 MHz with an electro-optic modulator (EOM). The beam reflected by the cavity is then detected by the photodiode Ph r . This signal is used to drive the acousto-optic modulator (AOM) frequency for a fast control and the Peltier element of the laser for a slow control of the laser frequency.

Our birefringence measurement is based on an ellipticity measurement. Light is polarized just before entering the cavity by the polarizer P. The beam transmitted by the cavity is then analyzed by the analyzer A crossed at maximum extinction and collected by a low noise photodiode Ph e (intensity of the extraordinary beam I e ). The analyzer also has an escape window which allows us to extract the ordinary beam (intensity I t ) which corresponds to the polarization parallel to P. This beam is collected by the photodiode Ph t .

All the optical components from the polarizer P to the analyzer A are placed in a ultra high vacuum chamber. In order to perform birefringence measurements on high purity gases, the vacuum chamber is connected to several gas bottles through leak valves which allow to precisely control the amount of injected gas. Finally, since the goal of the experiment is to measure magnetic birefringence, magnets surround the vacuum pipe. The transverse magnetic field is created thanks to pulsed coils described in Ref. [16] and briefly detailed in the next section.

Both signals collected by the photodiodes outside the cavity are simultaneously used in the data analysis as follows:

I e I t = σ 2 + Ψ 2 tot , (3)
where Ψ tot is the total ellipticity acquired by the beam going from P to A and σ 2 is the polarizer extinction ratio. Our polarizers are Glan laser prisms which have an extinction ratio of 2 × 10 -7 . The origin of the total ellipticity of the cavity is firstly due to the intrinsic birefringence of the mirrors M 1 and M 2 , as it will be discussed in section II C 2. We define the ellipticity imparted to the linearly polarized laser beam when light passes trough each mirror substrate as Γ s1,2 , and the one induced by the reflecting layers of the mirrors as Γ c . An additional component Ψ of the total ellipticity can be induced by the external magnetic field. Since we use pulsed magnetic fields, this ellipticity is a function of time.

Finally, if the ellipticities are small compared to unity, one gets:

I e (t) I t (t) = σ 2 + [Γ + Ψ(t)] 2 , (4)
where Γ = Γ s1 + Γ s2 + Γ c is the total static birefringence.

B. Magnetic field

It is clear from Eq. ( 1) that one of the critical parameter for experiments looking for magnetic birefringence is B 2 L. Our choice has been to reach a B 2 L as high as possible having a B as high as possible with a L such as to set-up a table-top low noise optical experiment. This is fulfilled using pulsed magnets that can provide fields of several tens of Tesla. Our apparatus consists of two magnets, called Xcoils. The principle of these magnets and their properties are described in details in Refs. [16,17].

The magnetic field profile along the longitudinal z-axis, which corresponds to the axis of propagation of the light beam, has been measured with a calibrated pick-up coil. Fig. 2 shows the normalized profile of an Xcoil. The magnetic field is not uniform along z. We define B max as the maximum field provided by the coil at its center and L B as the equivalent length of a magnet producing a uniform magnetic field B max such that:

+∞ -∞ B 2 (z)dz = B 2 max L B . (5)
L B is about the half of the Xcoil's length. Each Xcoil currently used has reached more than 14 T over 0.13 m of effective length corresponding to 25 T 2 m. The total duration of a pulse is a few milliseconds. The magnetic field reaches its maximum value within 2 ms. FIG. 2: Normalized profile of the square of the magnetic field along the longitudinal z-axis (solid line) inside one Xcoil. This is compared to the equivalent uniform magnetic field (dashed line) over the effective magnetic length LB (see text).

The pulsed coils are immersed in a liquid nitrogen cryostat to limit the consequences of heating which could be a cause of permanent damage of the coil's copper wire. Pulse duration is short enough that the coil, starting at liquid nitrogen temperature, always remains at a safe level i.e. below room temperature. A pause between two pulses is necessary to let the magnet cool down to the equilibrium temperature which is monitored via the Xcoils' resistance. The maximum repetition rate is 5 pulses per hour.

C. Fabry-Perot cavity

The other key point of our experiment is to accumulate the effect due to the magnetic field by trapping the light between two ultra high reflectivity mirrors constituting a Fabry-Perot cavity. Its length has to be large enough to leave a wide space so as to insert our two cylindrical cryostats (diameter of 60 cm for each cryostat) and vacuum pumping system. The length of the cavity is L c = 2.27 m corresponding to a free spectral range of F SR = c/2nL c ≃66 MHz with n the index of refraction of the considered medium in which the cavity is immersed. This index of refraction can be considered equal to one. The total acquired ellipticity Ψ is linked to the ellipticity ψ acquired in the absence of cavity, and depends on the cavity finesse F as follows [20]:

Ψ = 2F π ψ, (6)
where F is given by:

F = π √ R M 1 -R M , (7)
with R M the intensity reflection coefficient supposed to be the same for both mirrors. In order to increase the induced signal, a finesse as high as possible is essential.

Cavity finesse and transmission

Experimentally, the finesse is inferred from a measurement of the photon lifetime τ inside the cavity as presented on Fig 3. For t < t 0 , the laser is locked to the cavity. The laser intensity is then switched off at t 0 thanks to the AOM shown on Fig. 1 and used as an ultrafast commutator. For t > t 0 , one sees the typical exponential decay of the intensity of the transmitted ordinary beam [START_REF] Svelto | Principles of lasers[END_REF]:

I t (t) = I t (t 0 )e -(t-t0)/τ . (8)
The photon lifetime is related to the finesse of the cavity through the relation:

τ = nL c F πc . (9)
By fitting our data with Eq. ( 8) we get τ = 1.16 ms corresponding to a finesse of F = 481 000 and a cavity linewidth of ∆ν = c/2nL c F = 137 Hz. We summarize in Table I the performances of some well-known sharp cavities at λ = 1064 nm. To our knowledge, currently our interferometer is the sharpest in the world.

The transmission of the cavity T c is another important parameter. It corresponds to the intensity transmitted by the cavity divided by the intensity incident on the cavity when the laser frequency is locked. Indeed in order not I: Performances summary of the sharpest infra-red interferometers in the world. Lc is the length of the Fabry-Perot cavity, F SR is its full spectral range, F is the cavity finesse, τ is the photon lifetime, ∆ν is the frequency linewidth and Q = ν laser /∆ν is the quality factor of the interferometer, with ν laser the laser frequency. to be limited by the noise of photodiodes Ph t and Ph e , I t and I e have to be sufficiently high. This point is particularly critical for I e which corresponds to the intensity transmitted by the cavity multiplied by σ 2 . With a Ph e noise equivalent power of 11 fW/ √ Hz, we need an incident power greater than 0.2 nW so as not to be limited by electronic noise of Ph e .

Our cavity transmission is 20 %. The measurements of the finesse and the transmission allow to calculate mirrors properties such as their intensity transmission T M and their losses P M thanks to the following relations:

F = π T M + P M , (10)
T c = T M F π 2 , (11)
supposing that the mirrors are identical. We found T M = 3 ppm and P M = 3.5 ppm, which corresponds to the specifications provided by the manufacturer.

To conclude, our high finesse cavity enhances the Cotton-Mouton effect of a factor 2F/π=306 000, and its transmission allows measurements that are not limited by the noise of the detection photodiodes.

Cavity birefringence

The origin of the total static ellipticity is due to the mirror intrinsic phase retardation. Mirrors can be regarded as wave plates and for small birefringence, combination of both wave plates gives a single wave plate. The phase retardation and the axis orientation of this equivalent wave plate depend on the birefringence of each mirror and on their respective orientations [19,20].

The intrinsic phase retardation of the mirrors is a source of noise limiting the sensitivity of the apparatus. Moreover, since our signal detection corresponds to a homodyne technique, the static ellipticity Γ is used as a zero frequency carrier. To reach a shot noise limited sensitivity, one needs Γ to be as small as possible [17], implying that the phase retardation axis of both mirrors have to be aligned. For magnetic birefringence measurements, both mirrors' orientation is adjusted in order to have 10 -3 < Γ < 3 × 10 -3 rad.

Measurement of the total ellipticity as a function of mirror orientation allows to calculate the mirror intrinsic phase retardation per reflection. The experimental procedure is presented in Ref. [26]. The deduced phase retardation for our mirrors is δ M = (7 ± 6) × 10 -7 rad. Although the origin of the mirrors' static birefringence is still unknown, a review of existing data shows that for interferential mirrors phase retardation per reflection decreases when reflectivity increases [26]. This observation is confirmed by our new measurement. It is also in agrement with the empirical trend given in Ref. [26]:

δ M ≃ 0.1 × (1 -R M )
. Numerical calculations show that this trend can be explained assuming that the effect is essentially due to the layers close to the substrate.

As said before, mirror birefringence has two contributions: one comes from the substrate while the other is due to the reflecting layers. Whereas previous measurements do not allow to distinguish between the two contributions, we will see that this can be achieved with the measurement of I e decay.

A typical time evolution of I e when the incident beam locked to the cavity is switched off is shown in Fig. 4. We see that this curve can not be fitted by an exponential decay. As explained in Ref. [25], one has to take into account the intrinsic birefringence of the cavity. Nevertheless, the expression derived in Ref. [25], which only takes into account the reflecting layers birefringence, does not always fit our data. The evolution of I e presents sometimes an unexpected behavior: whereas no photon enters anymore into the cavity at t = t 0 , the extraordinary intensity starts growing before decreasing. To reproduce this behavior, one has to take into account the substrate birefringence. Lets calculate the transmitted intensity along the round-trip inside the cavity:

• For t ≤ t 0 , the laser is continuously locked to the cavity. According to Eq. ( 4), the intensities of the ordinary and the extraordinary beams are related by:

I e (t ≤ t 0 ) = [σ 2 + (Γ s2 + Γ s1 + Γ c ) 2 ]I t (t ≤ t 0 ). ( 12
)
• At t = t 0 , the laser beam is abruptly switched off, the cavity empties gradually. The ordinary and extraordinary beams are slightly transmitted at each reflection on the mirrors. But, because these mirrors are birefringent, some photons of the ordinary beam are converted into the extraordinary one. The reverse effect is neglected because I e ≪ I t .

We then follow the same procedure as in Ref. [25] to calculate the time evolution of I e . For t > t 0 , one gets:

I e (t) = ( 13
)
I t (t) σ 2 + Γ s1 + Γ s2 + Γ c 1 + t -t 0 2τ 2 .
The behavior shown on Fig.

4 is reproduced if (Γ s1 + Γ s2 ) ≃ -Γ c .
This expression is used to fit our experimental data plotted on Fig. 4. We find a photon lifetime of τ = 960 µs which is in good agreement when fitting I t [START_REF] Fig | run of data and thus can't be directly compared[END_REF], (Γ s1 + Γ s2 ) = 2 × 10 -3 rad and Γ c = -7 × 10 -3 rad.

We have for the first time the evidence that the substrate is birefringent and that this birefringence contributes to the total ellipticity due to the cavity.

D. Signal analysis

The starting point of our analysis are the voltage signals V e and V t provided by Ph e and Ph t . Voltage signals have to be converted into intensity signals by using the photodiode conversion factor g e and g t :

I e = g e V e , (14)
I t = g t V t . (15)
As demonstrated in Ref. [25], before analyzing raw signals one has to take into account the first order low pass filtering of the cavity. I t,filtered in the Fourier space is given by:

I t,filtered (ω) = 1 1 + i ω ωc I t (ω), (16)
where ν c = ω c /2π = 1/4πτ is the cavity cutoff frequency. Then, according to Eq. ( 4), the ellipticity Ψ(t) to be measured can be written as:

Ψ(t) = -Γ + I e (t) I t,filtered (t) -σ 2 . ( 17
)
The total static birefringence Γ is measured a few milliseconds just before the beginning of the magnetic pulse, thus when Ψ(t) = 0. On the other hand, Ψ is proportional to the square of the magnetic field and thus can be written as:

Ψ(t) = κB 2 filtered (t). ( 18
)
Since the photon lifetime is comparable with the rise time of the magnetic field, the first-order low pass filtering of the cavity has also to be taken into account on the quantity B 2 (t) as in Ref. [25]. To recover the value of the constant κ we calculate for each pulse the correlation between Ψ(t) and B 2 filtered (t):

κ = Ti 0 Ψ(t)B(t) 2 filtered dt Ti 0 [B(t) 2 filtered ] 2 dt , (19)
where T i is the integration time. A statistical analysis gives the mean value of κ and its uncertainty.

The magnetic birefringence ∆n is finally given by:

∆n(T, P ) = κ 4π τ F SR × λ L B × 1 sin 2θ . ( 20
)
∆n is thus expressed in T -2 . T and P correspond to gas temperature and pressure when measurements of magnetic birefringence on gases are performed. We define the normalized birefringence ∆n u as ∆n for P = 1 atm and B = 1 T.

III. EXPERIMENTAL PARAMETERS AND ERROR BUDGET

In the following, to evaluate the precision of our apparatus in the present version, we list the uncertainties at 1σ on the measurement of the parameters of Eq. ( 20) as recommended in Ref. [START_REF] Mohr | [END_REF]. The uncertainty on the magnetic birefringence has two origins. The evaluation of the uncertainty by a statistical analysis of series of observations is termed a type A evaluation and mainly concerns the measurement of τ and κ. An evaluation by means other than the statistical analysis of series of observations, calibrations for instance, is termed a type B evaluation and especially affects the parameters B, F SR, L B , λ and θ.

A. Photon lifetime in the Fabry-Perot Cavity

The photon lifetime τ is measured by analyzing the exponential decay of the intensity of the transmitted light. Several measurements have been performed both before and after almost each magnetic pulse. The uncertainty on the value of τ comes from the fact that mirrors can slightly move because of thermal fluctuations and acoustic vibrations. Measurements conducted in the same experimental conditions have been studied statistically leading to a relative variation of τ that does not exceed 2 % at 1σ-level. Data taken during operation, i.e. before and after magnetic pulses, show the same statistical properties as the ones taken without any magnetic field. Thus, the magnetic field does not cause additional change in τ .

B. Correlation factor

The correlation factor κ is given by Eq. ( 19). The A-type uncertainty on κ depends on the measurement of Ψ and thus on the experimental parameters given in Eq. (17). In practice, we pulse the magnets several times in the same experimental conditions to obtain a set of values of κ. The distribution of the κ values is found to be gaussian, and we assume that its standard deviation corresponds to the A-type uncertainty on κ. For our measurements performed with nitrogen and presented in section IV B, the A-type relative uncertainty is typically 3.5 %. The standard uncertainty of the average value of κ can then be reduced increasing the number of pulses.

B-type uncertainties depend on those of the square of the magnetic field, the photodiode conversion factors, and the filter function applied to the field.

To measure the magnetic field during operation, we measure the current which is injected in our X-coil. As mentioned in Ref. [16], the form factor B/I has been determined experimentally during the test phase by varying the current inside the X-coil (modulated at room temperature or pulsed at liquid nitrogen temperature), and by measuring the magnetic field induced on a calibrated pick-up coil. These measurements have led to a relative B-type uncertainty of δB/B = 0.7 % for the magnetic field corresponding to a B-type uncertainty on κ of 1.4 %.

The ratio g e /g t is measured from time to time by sending the same light intensities to each photodiode. The relative uncertainty in this parameter is 1.5 % corresponding to the same amount relative uncertainty in κ.

I t (t) and B 2 (t) are also filtered by a function that involves the parameter τ . We have empirically determined that a τ -variation of 2% led to a κ-variation of 0.8 %.

We can finally add quadratically the uncertainties above, and deduce that a B-type uncertainty of 2.2 % must be taken into account on every measurement of the correlation factor κ.

C. Frequency splitting between perpendicular polarizations

In this section we evaluate the attenuation of the extraordinary beam transmitted by our sharp resonant Fabry-Perot cavity on which the laser's ordinary beam is frequency-locked. Let's suppose that the ordinary (resp. extraordinary) beam resonates in the interferometer at the frequency ν t (resp. ν e ). The laser is locked to the cavity thanks to the ordinary beam. Thus ν t corresponds to the top of the transmission Airy function A of the Fabry-Perot cavity which is given by:

A(ν) = T c 1 + 4F 2 π 2 sin 2 2πnLc c ν . (21)
The frequency ν e is shifted from ν t by a quantity δν as it is shown on Fig. 5. The frequency splitting δν = ν t -ν e can be expressed as a function of the phase retardation δ acquired along a round-trip between the ordinary and the extraordinary beams:

δν = c 2πnL c δ, = F ∆ν π δ. ( 22
)
This formula indicates that in order to have a splitting very small compared to the cavity linewidth (δν ≪ ∆ν), the phase retardation δ must satisfy the following condition:

δ ≪ π F , (23)
which is equivalent to the condition on the acquired total ellipticity Ψ:

Ψ ≪ 1. ( 24
)
By combining Eqs. ( 21) and ( 22), we obtain the factor of attenuation a of the transmitted extraordinary beam's intensity given by:

a = A(ν e ) A(ν t ) , = 1 1 + 4F 2 π 2 sin 2 2πnLc c δν , = 1 1 + 4F 2 π 2 sin 2 δ . ( 25
)
The attenuation factor a is plotted as a function of δ on Fig. 6 The frequency splitting can first be due to our birefringent cavity. As in Ref. [20], let's consider both cavity mirrors equivalent to a single wave-plate with phase retardation δ w = δ between both polarizations. The total phase retardation δ w is linked to the cavity mirrors' M 1 and M 2 own phase retardation δ 1 and δ 2 as [20]:

δ w = (δ 1 -δ 2 ) 2 + 4δ 1 δ 2 cos 2 (θ m ). ( 26
)
To set a δ w as small as possible so as to minimize the correction to I meas e , one needs to adjust the angle θ m between the neutral axes of both mirrors. This way, we set a δ w of the order of a few 10 -8 rad, corresponding to a correction smaller than 0.001 % on I meas e . Secondly, the frequency splitting between both polarizations can be due to the induced magnetic birefringence of the medium inside the chamber. As seen above, the induced ellipticity given by Eq. ( 24) must be well below 1 rad. This condition is always satisfied in the range of pressure and field we are working. The induced ellipticity does not exceed 10 -2 rad. This corresponds at worst to a phase retardation of δ = 10 -7 rad. The attenuation factor I meas e is thus smaller than 0.1 %. In principle, this attenuation generates an error that has to be taken into account in the measured ratio I e /I t,filtered of Eq. ( 17), which implies an error in the value of κ. At present, since the attenuation is smaller than 0.1 %, this error can be neglected compared to the others uncertainties in κ.

D. Cavity free spectral range

The dedicated experimental setup for the measurement of the cavity free spectral range F SR = c/2nL c is shown on Fig. 7. The principle is to inject into the cavity two laser beams shifted one compared with the other by a given frequency. This frequency is then adjusted to coincide with the free spectral range.

Experimentally, the main beam is divided into two parts thanks to a polarizing beam splitting cube. The first part is directly injected into the cavity while the other one is frequency shifted by the acousto-optic modulator AOM2 with a double-pass configuration before injection. The main beam is frequency modulated with a voltage ramp applied on a piezo element mounted on the crystal resonator of the laser. The intensity transmitted by the cavity is observed on Ph t as shown on Fig. 8. The solid line corresponds to the intensity of the first beam. We observe typical Fabry-Perot peaks whose frequency gap corresponds to F SR. Peaks due to the second beam (dashed line) are frequency shifted by 2f AOM2 . We finally adjust f AOM2 in order to superimpose both series of peaks. The precise knowledge of the driven frequency f AOM2 enables us to determine with the same precision the value of the free spectral range, and thus the cavity length. The adjustment of fAOM2 in order to superimpose both series of peaks allows to precisely measure the free spectral range F SR of the cavity.

A typical value is F SR = (65.996 ± 0.017) MHz. This corresponds to a cavity length of L c = (2.2713 ± 0.0006) m. Since this length can be prone to variation, the F SR value is regularly checked and updated.

E. Effective magnetic length

Following Eq. ( 5), the effective magnetic length L B has been calculated by numerically integrating the field measured with a calibrated pick-up coil. Taking into account the experimental uncertainties, we got for one Xcoil: L B = (0.137 ± 0.003) m, corresponding to a relative B-type uncertainty on L B of 2.2 %.

F. Laser wavelength

As mentioned above, infra-red light enters the cavity. The wavelength of the Nd:YAG laser is 1064 nm, and its uncertainty is given by the width of the laser transition. The natural linewidth of Nd:YAG lasers are not usually given by the manufacturers. However, we can estimate it from the bandwidth of the gain curve of the amplifying medium. It is typically of the order of 30 GHz [START_REF] Hecht | The Laser Guidebook[END_REF]. This corresponds to an uncertainty on the laser wavelength of 0.3 nm. In order to be conservative, we use λ = (1064.0± 0.5) nm. The relative uncertainty is negligible in our case, compared to main uncertainties.

G. Angle between the incident polarization and the magnetic field direction

The angle between the incident light polarization and the magnetic field direction is adjusted to 45 • thanks to magnetic birefringence measurements as a function of the polarizer direction θ P . In order to be more sensitive, this is performed close to the position where the magnetic field is parallel to the polarizer P (θ = 0 • ).

Measurements are realized with about 7 × 10 -3 atm of air. The analyzer direction is crossed at maximum extinction each time the polarizer is turned. Fig. 9 represents the evolution of the correlation factor κ as a function of θ P . Data are fitted by a sinusoidal trend κ(θ P ) = κ 0 sin 2(θ P -θ 0 ) giving θ 0 = (2.6 ± 0.2) • . This measurement allows to set θ to (45.0 ± 1.2) • . The uncertainty is mainly due to the mechanical system which holds and turns the polarizer.

H. Error budget

We summarize in the Table II the typical values of the experimental parameters that have to be measured and their B-type associated uncertainty. These uncertainties are quadratically added to give a B-type relative uncertainty on the birefringence ∆n of 3.1 % at 1σ. I. Temperature and pressure of gases magnetic birefringence measurements are performed at room temperature T = 293 K. The experimental room is air-conditioned. A flow of compressed air between the outer wall of the vacuum pipe and the liquid nitrogen cryostat containing the magnet maintains the room temperature in the gas chamber.

A temperature profile has been realized along the length of the vacuum pipe, and is plotted on Fig. 10. The temperature variation does not exceed 1 K inside the tube that passes through the magnetic field. Concerning gases, we consider that our birefringence measurements are given at (293 ± 1) K. The pressure of the gas inside the chamber is measured at each side of the vacuum pipe getting into magnets with pressure gauges. The relative uncertainty provided by the manufacturer is 0.2 %.

IV. MAGNETIC BIREFRINGENCE MEASUREMENTS

A. Raw signals The laser beam remains locked to the Fabry-Perot cavity, despite mechanical vibrations caused by the shot of magnetic field. I t and I e start oscillating after about 4 ms. Seismometers placed on mirror mounts show that these oscillations are mainly due to acoustic perturbations produced by the magnet pulse and propagating from the magnet to mirror mounts through the air. We also see that the minimum of I e does not coincide with the maximum of B 2 . This phenomenon is due to the cavity filtering as explained in details in Ref. [25]. In Fig. 12, we plot the square of the magnetic field filtered by the cavity and the ellipticity calculated with Eq. ( 17) as a function of time. Since the acoustic perturbations affect both signals I t and I e , and taking into account the cavity filtering between I t and I e , oscillations on Ψ are strongly reduced to a few 10 -5 rad, thus not visible on this figure. These oscillations induce uncertainty to the measurement but are already included in the Atype uncertainty on kappa measured in section III B.

Finally, we note that both quantities B 2 filtered and Ψ reach their extremum at the same time and their variation can be perfectly superimposed, providing a very precise measurement of magnetic linear birefringence of nitrogen gas. FIG. 12: Cotton-Mouton effect measurement on 32.1 × 10 -3 atm of molecular nitrogen. Gray line: Total ellipticity as a function of time. Dashed line: Square of the magnetic field filtered by a first-order low pass filter corresponding to the cavity filtering.

B. Apparatus calibration

In order to calibrate our apparatus and to evaluate its present sensitivity we have measured the magnetic birefringence of molecular nitrogen. These measurements have been performed at different pressure from 2.1 × 10 -3 to 32.1 × 10 -3 atm and are summarized in Fig. 13. In this range, nitrogen can be considered as an ideal gas and the pressure dependence of its birefringence is thus linear:

∆n(T -2 ) = ∆n u (atm -1 T -2 ) × P (atm). ( 27
)
We have checked that our data are correctly fitted by a linear equation. Its ∆n axis-intercept is consistent with zero within the uncertainties. Its slope gives the normalized magnetic birefringence at B = 1 T and P = 1 atm:

∆n u = (-2.00 ± 0.08 ± 0.06) × 10 -13 atm -1 T -2 .
The first uncertainty 0.08 × 10 -13 atm -1 T -2 corresponds to the fitting uncertainty and represents the Atype total uncertainty at 1σ. The second one 0.06 × 10 -13 atm -1 T -2 represents the B-type uncertainty at 1σ.

Our value of the normalized birefringence is compared in Table . III to other experimental published values at λ = 1064 nm [31,32]. This shows that our value perfectly well with other existing measurements. Our total uncertainty is 10 -14 atm -1 T -2 , calculated by quadratically adding the A-type and B-type uncertainties. This is 1.8 times more precise than the other results. It therefore provides a successful calibration of the whole apparatus.

Ref.

∆nu × 10 -13 (at P = 1atm and B = 1T) [31] -2.17 ± 0.21 [32] -2.02 ± 0.

C. Upper limit on vacuum magnetic birefringence measurements

Once the calibration performed we have evaluated the upper limit of the present apparatus on vacuum magnetic birefringence. To this end, several pulses were performed in vacuum. In Fig. 14, a typical ellipticity measured during a magnetic pulse is plotted. Acoustic perturbations induce oscillations of Ψ starting at about 4 ms, with variations of the order of 10 -5 rad. In order to infer our best upper limit for the value of the vacuum magnetic birefringence, we limit the integration time to 4 ms. We get ∆n < 5.0 × 10 -20 T -2 per pulse. During operation, the pressure inside the UHV system was better than 10 -10 atm. To be conservative, lets residual gases are mainly 78 % of nitrogen and 21 % of oxygen. The normalized magnetic birefringences of these gases are of the order of -2 × 10 -13 atm -1 T -2 and -2 × 10 -12 atm -1 T -2 respectively [6]. The total residual magnetic birefringence is then of the order of 6 × 10 -23 T -2 , which is well below our current upper limit. In the final setup, vacuum quality will be monitored with a residual gas analyzer.

V. CONCLUSION

The successful calibration we report in this paper is a crucial step towards the measurement of the vacuum magnetic birefringence. It shows our capability to couple intense magnetic fields with one of the sharpest Fabry-Perot cavity in the world. It is worthwhile to note that an energy of about 100 kJ is discharged in our coils during a few milliseconds. These 10 MW of electrical power generate acoustic perturbations and mechanical vibrations that tend to misalign the cavity mirrors. The linewidth ∆ν of our Fabry-Perot cavity is of the order of 150 Hz. A relative displacement ∆L c = L c × ∆ν/ν laser = 1 pm of both mirrors is enough to get out of resonance. The sharper the cavity, the bigger the challenge.

The sensitivity per pulse we got both in gases and in vacuum is the best ever reached for this kind of measurement. For sake of comparison, the best birefringence limit obtained in vacuum with continuous magnets is ∆n ≤ 2.1 × 10 -20 T -2 with an integration time of t int = 65 200 s [13]. In order to compare both methods, we need to translate the best limit obtained in continuous regime to the one obtained with our integration time T i = 4 ms. Assuming white noise for both methods, best limit reported in Ref. [13] corresponds to ∆n(T i ) = ∆n(t int ) t int /T i ≤ 8.5 × 10 -17 in 4 ms of integration. This value is more than three orders of magnitude higher than ours, proving that pulsed fields are a powerful tool for magnetic birefringence measurements.

Long terms perspective is to get a value of ∆n = 4 × 10 -24 T -2 , corresponding to the vacuum magnetic birefringence, with at most 1000 pulses. This corresponds to a sensitivity better than 1.3 × 10 -22 T -2 per pulse. A factor of the order of 10 on optical sensitivity will be achievable with a better acoustic insulation and a more robust locking system, especially reducing the noise of the measured light intensities transmitted by the cavity. Further improvements depend on the possibility to have higher magnetic fields. We have designed a new pulsed coil, called XXL-coil, which has already reached a field higher than 30 T when a current higher than 27 000 A is injected. This corresponds to more than 300 T 2 m [33]. Two XXL-coils will allow us to improve our current sensitivity by a factor 100. In the near future, the apparatus will be modified in order to host these XXL-coils. Therefore the final version of the experiment will be ready for operation.

4 BFIG. 1 :

 41 FIG.1: Experimental setup. A Nd-YAG laser is frequency locked to the Fabry-Perot cavity consisting of the mirrors M1 and M2. The laser beam is linearly polarized by the polarizer P and analyzed with the polarizer A. This analyzer allows to extract the extraordinary beam sent on photodiode Phe as well as the ordinary beam sent on photodiode Pht. The beam reflected by the cavity analyzed on the photodiode Phr is used for the cavity locking. A transverse magnetic field B can be applied inside the cavity in order to study the magnetic birefringence of the medium. EOM = electro-optic modulator; AOM = acousto-optic modulator, PDH = Pound-Drever-Hall.

FIG. 3 :

 3 FIG.3: Time evolution of the intensity of the ordinary beam (gray line). The laser is switched off at t = t0. Experimental data are fitted by an exponential decay (black dashed line) giving a photon lifetime of τ = 1.16 ms, a finesse of F = 481 000 and a linewidth of ∆ν = c/2nLcF = 137 Hz.

FIG. 4 :

 4 FIG.4: Time evolution of the intensity of the extraordinary beam (gray line). The laser is switched off at t = t0. Experimental data are perfectly fitted by Eq. (13) (black dashed line).

FIG. 5 :

 5 FIG.5: Airy function of our Fabry-Perot cavity (linewidth ∆ν = 137 Hz and transmission Tc = 20 %). The frequency of the ordinary beam is assumed to be locked at the top of the transmission function (solid line) while he frequency νe of the extraordinary beam is shifted from νt of a quantity δν (dashed line).

  for a finesse F = 481 000. The real intensity I e of the extraordinary beam transmitted by the cavity is obtained from the corrected measured intensity I meas e as I e = I meas e / a.

FIG. 6 :

 6 FIG.6: Attenuation factor a as a function of the phase retardation δ between both polarizations.

FIG. 7 :

 7 FIG.7: Experimental setup for the cavity length measurement. Two laser beams frequency shifted one compared to the other by AOM2 are injected into the cavity. The frequency of the laser is frequency modulated with a voltage ramp applied on a piezo element mounted on the crystal resonator of the laser. Photodiode Pht allows to observe the typical Fabry-Perot peaks from which the F SR measurement is performed.

FIG. 8 :

 8 FIG.8: Transmission peaks of the Fabry-Perot cavity as a function of the laser frequency. Two beams are sent to the interferometer: the second beam (dashed line) is frequency shifted by 2fAOM2 compared to the first beam (solid line). The adjustment of fAOM2 in order to superimpose both series of peaks allows to precisely measure the free spectral range F SR of the cavity.

FIG. 9 :

 9 FIG.9: Correlation factor κ between the square of the magnetic field and the ellipticity as a function of the angle θP of the incident polarization.

FIG. 10 :

 10 FIG.10: Profile of the temperature inside the vacuum pipe along the longitudinal z-axis. The X-coil is also schematized at the center. The temperature variation does not exceed 1 K inside the tube that runs through it.

Fig. 11 presents

 11 Fig. 11 presents signals obtained with 32.1 × 10 -3 atm of molecular nitrogen. The intensity of the ordinary beam I t (top) remains almost constant while the intensity of the extraordinary beam I e (middle) varies when the magnetic field (bottom) is applied. The magnetic field reaches its maximum of 5.2 T within less than 2 ms.The laser beam remains locked to the Fabry-Perot cavity, despite mechanical vibrations caused by the shot of magnetic field. I t and I e start oscillating after about 4 ms. Seismometers placed on mirror mounts show that these oscillations are mainly due to acoustic perturbations produced by the magnet pulse and propagating

  FIG. 11:Cotton-Mouton effect measurement on 32.1 × 10 -3 atm of molecular nitrogen.(top) Normalized intensity of the ordinary beam as a function of time. (middle) Intensity of the extraordinary beam divided by the mean of It as a function of time. (bottom) Square of the magnetic field as a function of time.

FIG. 13 :

 13 FIG.13: Magnetic birefringence of molecular nitrogen as a function of pressure. The solid line corresponds to the linear fit of experimental data.

FIG. 14 :

 14 FIG.14: Typical ellipticity (gray) measured during a magnetic pulse (black) performed in vacuum. Acoustic perturbations generate ellipticity oscillations starting at 4 ms.

TABLE

TABLE II :

 II Parameters that have to be measured to infer the value of the birefringence ∆n and their respective relative Btype uncertainty at 1σ.

Acknowledgments

We thank all the members of the BMV collaboration, and in particular J. Béard, J. Billette, P. Frings, J. Mauchain, M. Nardone, L. Recoules and G. Rikken for strong support. We are also indebted to the whole technical staff of LNCMI. We acknowledge the support of the Fondation pour la recherche IXCORE and of the ANR-Programme non thématique (ANR-BLAN06-3-139634).


Pierrick Coupé

Pierre Hellier

Xavier Morandi

Christian Barillot

3D Rigid Registration of Intraoperative Ultrasound and Preoperative MR Brain Images based on Hyperechogenic Structures

The registration of intraoperative ultrasound (US) images with preoperative magnetic resonance (MR) images is a challenging problem due to the difference of information contained in each image modality. To overcome this difficulty, we introduce a new probabilistic function based on the matching of cerebral hyperechogenic structures. In brain imaging, these structures are the liquid interfaces such as the cerebral falx and the sulci, and the lesions when the corresponding tissue is hyperechogenic. The registration procedure is achieved by maximizing the joint probability for a voxel to be included in hyperechogenic structures in both modalities. Experiments were carried out on real datasets acquired during neurosurgical procedures. The proposed validation framework is based on (i) visual assessment, (ii) manual expert estimations , and (iii) a robustness study. Results show that the proposed method (i) is visually efficient, (ii) produces no statistically different registration accuracy compared to manual-based expert registration, and (iii) converges robustly. Finally, the computation time required by our method is compatible with intraoperative use.

Introduction

Due to its low cost, its real-time imaging capabilities, and its non invasive nature, ultrasound (US) imaging has become a popular modality. These attributes have been used for a large number of clinical applications. In neurosurgery, 1 ultrasound imaging has been employed in many cases of brain examinations over the last two decades [START_REF] Rubin | Intraoperative ultrasound examination of the brain[END_REF]. Several studies demonstrated that ultrasonography can be used for locating tumors, defining their margins, differentiating their internal characteristics, and for detecting of brain shift and residual tumoral tissues [START_REF] Dohrmann | History of intraoperative ultrasound in neurosurgery[END_REF]. At present, 3D US imaging is integrated within the neuronavigation systems to provide a useful and efficient intraoperative tool [START_REF] Unsgaard | Intra-operative 3D ultrasound in neurosurgery[END_REF]. Ultrasound imaging has also been shown to be a promising method for quantifying and for correcting brainshift in Image-Guided Neurosurgery (IGNS) [START_REF] Arbel | Automatic Nonlinear MRI-Ultrasound Registration for the Correction of Intra-operative Brain Deformations[END_REF][START_REF] Bucholz | The correction of stereotactic inacuracy caused by brain shift using an intraoperative ultrasound device[END_REF][START_REF] Comeau | Intraoperative ultrasound for guidance and tissue shift correction in image-guided neurosurgery[END_REF][START_REF] Gobbi | Ultrasound/MRI Overlay with Image Warping for Neurosurgery[END_REF][START_REF] Hata | Registration of ultrasound echography for intraoperative use: a newly developed multiproperty method[END_REF][START_REF] King | Bayesian estimation of intra-operative deformation for image-guided surgery using 3-D ultrasound[END_REF][START_REF] Letteboer | Brain Shift Estimation in Image-Guided Neurosurgery Using 3D Ultrasound[END_REF][START_REF] Lunn | Displacement Estimation with Co-Registered Ultrasound for Image Guided Neurosurgery: A Quantitative In Vivo Porcine[END_REF][START_REF] Pennec | Tracking brain deformations in timesequences of 3D US images[END_REF][START_REF] Reinertsen | Validation of vessel-based registration for correction of brain shift[END_REF][START_REF] Reinertsen | Clinical validation of vessel-based registration for correction of brain-shift[END_REF].

During a neurosurgical procedure, the ultrasound probe is tracked by the neuronavigation system which computes the 3D positions and orientations of the B-scans. Matching between the intraoperative US images and the preoperative MR image is ensured by a rigid registration. In phantom [START_REF] Comeau | Intraoperative ultrasound for guidance and tissue shift correction in image-guided neurosurgery[END_REF] and animal studies [START_REF] Lunn | Displacement Estimation with Co-Registered Ultrasound for Image Guided Neurosurgery: A Quantitative In Vivo Porcine[END_REF], the matching accuracy between intraoperative B-scans and preoperative images has been quantified between 1.5mm and 3mm. Nevertheless, in clinical context, the matching error can reach 10mm (see Table 1). This error includes tool calibration errors (the position localizer and the US probe), tool localization errors (tracking system error), and registration errors from the neuronavigation system.

Registration approaches based on classical image similarity measures such as the Sum Square Difference (SSD), Mutual Information (MI) or Correlation Ratio (CR) are known to fail to robustly register MR and US images [START_REF] Roche | Rigid registration of 3D ultrasound with MR images: a new approach combining intensity and gradient information[END_REF]. Therefore, other options have been investigated.

• a) Landmark-based registration represents the majority of the approaches [START_REF] Comeau | Intraoperative ultrasound for guidance and tissue shift correction in image-guided neurosurgery[END_REF][START_REF] Gobbi | Ultrasound/MRI Overlay with Image Warping for Neurosurgery[END_REF][START_REF] King | Bayesian estimation of intra-operative deformation for image-guided surgery using 3-D ultrasound[END_REF][START_REF] Porter | Three-dimensional registration and fusion of ultrasound and MRI using major vessels as fiducial markers[END_REF][START_REF] Reinertsen | Validation of vessel-based registration for correction of brain shift[END_REF][START_REF] Reinertsen | Clinical validation of vessel-based registration for correction of brain-shift[END_REF]. The motivation is bound to the difficulty of finding a function matching US image intensities with MR image intensities. These methods are based on the matching of manually defined points [START_REF] Gobbi | Ultrasound/MRI Overlay with Image Warping for Neurosurgery[END_REF], lines representing the vascular system [START_REF] Comeau | Intraoperative ultrasound for guidance and tissue shift correction in image-guided neurosurgery[END_REF][START_REF] Porter | Three-dimensional registration and fusion of ultrasound and MRI using major vessels as fiducial markers[END_REF][START_REF] Reinertsen | Validation of vessel-based registration for correction of brain shift[END_REF][START_REF] Reinertsen | Clinical validation of vessel-based registration for correction of brain-shift[END_REF], or cortical surface [START_REF] King | Bayesian estimation of intra-operative deformation for image-guided surgery using 3-D ultrasound[END_REF].

• b) Intensity-based approaches using histogram-based similarity measures tend to overcome the problem by preprocessing the images in order to register more similar images [START_REF] Arbel | Automatic Nonlinear MRI-Ultrasound Registration for the Correction of Intra-operative Brain Deformations[END_REF][START_REF] Letteboer | Rigid registration of 3D ultrasound data of brain tumours[END_REF].

• c) By introducing the Bivariate Correlation Ratio (BCR), Roche et al. [START_REF] Roche | Rigid registration of 3D ultrasound with MR images: a new approach combining intensity and gradient information[END_REF] incorporated the transformation of MR to pseudo-US image as a function into the similarity measure.

In this paper, we propose a new objective function based on the matching of the cerebral hyperechogenic structures such as sulci and the cerebral falx, and the lesion when the corresponding tissue is hyperechogenic. The registration is achevied by maximizing the correlation value between the US image and the probabilistic map of hyperechogenic structures estimated from MR image. The proposed method is thus a compromise between landmark and intensity-based approaches:

• As with landmark-based approaches, only regions considered as relevant are used to drive the registration procedure. In our method, these regions are the hyperechogenic structures of the brain.

• As with intensity-based methods, the proposed approach does not require segmentation of the US image what is a challenging problem.

2 Materials and Methods

Method overview

The scheme of the overall workflow is presented in figure 1. First, the "hyperechogenic" structures present in MR image (i.e. the structures visible in MR image expected to be hyperechogenic in intraoperative US) are detected with the M Lvv operator [START_REF] Luc | Scale and the differential structure of images[END_REF][START_REF] Maintz | Evaluation of ridge seeking operators for multimodality medical image matching[END_REF]. In brain imaging, these structures are the liquid interfaces such as the cerebral falx and the sulci, in addition to the lesions when the corresponding tissue is hyperechogenic (e.g. carvernoma or glioma). The curvature-based M Lvv operator was first introduced in [START_REF] Luc | Scale and the differential structure of images[END_REF][START_REF] Maintz | Evaluation of ridge seeking operators for multimodality medical image matching[END_REF] before being used to detect the sulci and the cerebral falx in [START_REF] Goualher | Statistical sulcal shape comparisons: Application to the detection of genetic encoding of the central sulcus shape[END_REF][START_REF] Goualher | Three-dimensional segmentation and representation of cortical sulci using active ribbons[END_REF][START_REF] Goualher | Automated extraction and variability analysis of sulcal neuroanatomy[END_REF]. The US image and the probability map of the hyperechogenic structures extracted from MR image are then registered by maximizing the probability for a voxel to be included in hyperechogenic structures in both modalities. Contrary to histogram-based approaches that match all the information in both images, the proposed approach consists of matching only hyperechogenic structures [START_REF] Coupé | A probabilistic objective function for 3D rigid registration of intraoperative US and preoperative MR brain images[END_REF], which makes it more robust to artefacts such as acoustic shadows. Indeed, in US imaging, the bright areas provide information on the underlying structures whereas the dark areas can correspond to the underlying anatomical structure or acoustic shadows [START_REF] Hellier | Acoustic shadows detection, application to accurate reconstruction of 3D intraoperative ultrasound[END_REF]. Moreover, the accuracy of sulci matching is an important issue since these structures are used by the neurosurgeon during the neurosurgical procedure [START_REF] Pierre Jannin | Integration of sulcal and functional information for multimodal neuronavigation[END_REF]. Finally, by using the natural property of US imaging to detect the hyperechogenic structures, the method does not require segmentation of the US image. This way, the method is less sensitive to error of US image segmentation and is less time consuming during the intraoperative stage.

Probabilistic objective function

The proposed registration process is based on the estimation of the transformation T maximizing the joint probability for a voxel X = (x, y, z) to be included in hyperechogenic structures in both modalities:

T = arg max T Ω p(X ∈ Φ U S , T (X) ∈ Φ M R ) dX (1)
where p(X ∈ Φ U S ) is the probability for X to be included in an hyperechogenic structure from the US image and p(X ∈ Φ M R ) is the probability for X to be included in an hyperechogenic structure (in the sense of the ultrasound image) from the T1-w MR image. Assuming that the probabilities are independent, we can write:

T = arg max T Ω p(X ∈ Φ U S ).p(T (X) ∈ Φ M R ) dX (2)
Our objective function can be viewed as the maximization of the correlation value between the two probability maps of hyperechogenic structures extracted from both modalities.

Construction of the probability maps

In order to construct the probability maps, we define a function f matching the intensity of both the US image and the MR image with the probability for X to be included in hyperechogenic structures:

p(X ∈ Φ) = f (u(X)) (3)
where u : Ω -→ R is an image defined on Ω.

Intraoperative US image

For the intraoperative US image U , by definition f is the identity function.

p(X ∈ Φ U S ) = U (X) (4)
The intensity of U are only scaled between 0 and 1 during surgery to fit with our probabilistic framework.

Preoperative MR image

For the preoperative MR image V , the evaluation of f is done prior to surgery and is based on both the detection of the liquid interfaces with the M Lvv operator and the segmentation of the pathological tissues. The Lvv operator is a robust intensity-based curvature detector [START_REF] Luc | Scale and the differential structure of images[END_REF] based on the first and second derivatives of the image intensities. The first and second derivatives are combined to obtain an operator less sensitive to flat areas with low gradients. This kind of operator is used to detect ridge-like features in images, with negative value for crests in the intensity domain and positive value for valleys in the intensity domain. In [START_REF] Maintz | Evaluation of ridge seeking operators for multimodality medical image matching[END_REF], the M Lvv has been proposed for multimodal registration of CT and MR images. In our case, as in [START_REF] Goualher | Statistical sulcal shape comparisons: Application to the detection of genetic encoding of the central sulcus shape[END_REF][START_REF] Goualher | Three-dimensional segmentation and representation of cortical sulci using active ribbons[END_REF][START_REF] Goualher | Automated extraction and variability analysis of sulcal neuroanatomy[END_REF], the M Lvv is used to extract the hyperechogenic structures (sulci, cerebral falx) from T1-w MR image. Overall, curvature information has been used by several other authors to characterize cortical features [START_REF] Cachia | A generic framework for parcellation of the cortical surface into gyri using geodesic Voronoi diagrams[END_REF][START_REF] Rettmann | Automated sulcal segmentation using watersheds on the cortical surface[END_REF][START_REF] Tosun | Cortical reconstruction using implicit surface evolution: Accuracy and precision analysis[END_REF]. Most of these methods are based on geodesic curvature computed on cortical surfaces.

In T1-w MR images, the sulci are valleys (negative ridges) in the intensity domain. By using the positive values of M Lvv, the sulci and the cerebral falx can be efficiently detected [START_REF] Goualher | Statistical sulcal shape comparisons: Application to the detection of genetic encoding of the central sulcus shape[END_REF][START_REF] Goualher | Three-dimensional segmentation and representation of cortical sulci using active ribbons[END_REF][START_REF] Goualher | Automated extraction and variability analysis of sulcal neuroanatomy[END_REF]. Figures 2, 3 and4 show the positive values of M Lvv operator. Finally, our function f is defined as :

p(X ∈ Φ M R ) = M Lvv(V (X))I M 1 (X) + Ψ(X)I M 2 (X) (5)
where I M is the indicator function for the set M:

• M 1 = {X ∈ Ω, such that M Lvv(V (X)) 0}
• M 2 = {X ∈ Ω, such that X belongs to the lesional tissue}

As for the US intensities, the positive values of the M Lvv are scaled between 0 and 1. The M Lvv operator is defined in 3D as :

M Lvv(V (x, y, z)) = -1 2 w 2 ∂V (X) ∂x 2 ∂ 2 V (X) ∂y 2 + ∂ 2 V (X) ∂z 2 -2 ∂V (X) ∂y ∂V (X) ∂z ∂ 2 V (X) ∂y∂z + ∂V (X) ∂y 2 ∂ 2 V (X) ∂x 2 + ∂ 2 V (X) ∂z 2 -2 ∂V (X) ∂x ∂V (X) ∂z ∂ 2 V (X) ∂x∂z + ∂V (X) ∂z 2 ∂ 2 V (X) ∂x 2 + ∂ 2 V (X) ∂y 2 -2 ∂V (X) ∂x ∂V (X) ∂y ∂ 2 V (X) ∂x∂y (6)
where

w 2 = ∂V (X) ∂x 2 + ∂V (X) ∂y 2 + ∂V (X) ∂z 2
. Ψ(X) is the probability given to X in the segmentation of pathological tissue M 2 . Ψ is used to incorporate a priori on pathology. For pathological tissue such as cavernoma or low-grade glioma, Ψ(X) is high since these tissues are hyperechogenic.

Preprocessing of the MR data before surgery

First, skull stripping is performed from the T1-w MRI sequence [START_REF] Mangin | Robust brain segmentation using histogram scale-space analysis and mathematical morphology[END_REF]. We choose to remove the skull prior to M Lvv computation because this structure does not appear in the area of the craniotomy. The raw MR images are then denoised using an optimized Non Local Means filter1 [START_REF] Coupé | An Optimized Blockwise NonLocal Means Denoising Filter for 3-D Magnetic Resonance Images[END_REF] before applying the M Lvv operator to the brain tissues. The use of a denoising stage makes the computation of the M Lvv more stable. Indeed, the presence of noise may create false positive or negative curvatures which could bias the registration framework. After applying the M Lvv operator, only the positive values (i.e. the sulci and the falx) are kept in the processing stream. Finally, the M Lvv map and the segmentation M 2 are merged together (see Fig. 2, see Fig 3 and Fig. 4). In our experiments, the segmentation of pathology was manually performed by the neuroanatomist before the surgical procedure (see Fig. 1). The computational time required by preprocessing steps performed during preoperative stage was 4 minutes for skull stripping and 3 minutes for denoising on a Pentium M 2GHz. In addition, 3-8 minutes were required for manual segmentation of the lesion according to its size on a Stealth Station TREON (Medtronic Inc., Minneapolis, USA). Since these steps are performed before surgery, there is no impact of on practical value of the proposed method.

Data acquisition

T1-w SENSE 3D sequences were used to acquire preoperative T1-weighted MR images on a 3T Philips Gyroscan scanner (Best, the Netherlands). During the neurosurgical procedure, the US probe (Sonosite Inc. Bothell, WA. USA, cranial 7 -4M Hz probe) was tracked by the Polaris cameras of the Stealth Station TREON (Medtronic Inc., Minneapolis, USA). The Sononav software of the neuronavigation system was used to acquire the 2D B-scans and the probe positions. From the 2D B-scans and their positions, a 3D volume was reconstructed with the Probe Trajectory method [START_REF] Coupé | Probe Trajectory Interpolation for 3D Reconstruction of Freehand Ultrasound[END_REF]. The experiments were carried out on 3 patients. For each patient, a sequence of images was acquired before opening the dura. Some studies have considered quantitative measurement of brain shift during surgical procedures and showed that non significant displacement occurred before dura opening [START_REF] Hill | Measurement of introperative brain surface deformation under a craniotomy[END_REF][START_REF] Roberts | Intraoperative brain shift and deformation: a quantitative analysis of cortical displacement in 28 cases[END_REF]. Thus, we assumed that the transformation between intraoperative US and preoperative MR was rigid. The characteristics of reconstructed volumes are:

• for patient 1 a 3D volume of 486 × 462 × 206 voxels with a resolution of 0.15 × 0.14 × 0.14 mm 3 ,

• for patient 2 a 3D volume of 510 × 423 × 174 voxels with a resolution of 0.21 × 0.19 × 0.20 mm 3 ,

• for patient 3 a 3D volume of 265 × 450 × 324 voxels with a resolution of 0.19 × 0.17 × 0.18 mm 3 .

MR-US registration of the neuronavigation system

During all the neurosurgical procedure, the coordinate system of the preoperative MR image and the coordinate system of the intraoperative field are related by a rigid registration. The rigid registration of the neuronavigation system is based on surface matching between the preoperative MR image and the position of points acquired on the patient's head with the position localizer. First, the skin is extracted from the MR image by manual thresholding. A cloud of points is then continuously acquired on the patient's head close to the eyes region by moving the position localizer. Following this, one point is acquired on each ear with another point on the extremity of the patient's nose. Finally, the neuronavigation system performs a points to surface matching. According to phantom and animal studies, the errors in probe calibration, 3D localization of the probe, and rigid registration performed by the neuronavigation system lead to a global error less than 3mm [START_REF] Comeau | Intraoperative ultrasound for guidance and tissue shift correction in image-guided neurosurgery[END_REF][START_REF] Letteboer | Brain Shift Estimation in Image-Guided Neurosurgery Using 3D Ultrasound[END_REF][START_REF] Lunn | Displacement Estimation with Co-Registered Ultrasound for Image Guided Neurosurgery: A Quantitative In Vivo Porcine[END_REF]. The error due to the 3D localization of the probe is estimated to 0.35mm for each marker on a tool from the manufacturer [START_REF]northern digital inc[END_REF]. The error due to the calibration is generally estimated around 1.5mm [START_REF] Comeau | Intraoperative ultrasound for guidance and tissue shift correction in image-guided neurosurgery[END_REF][START_REF] Lunn | Displacement Estimation with Co-Registered Ultrasound for Image Guided Neurosurgery: A Quantitative In Vivo Porcine[END_REF]. In our case, the probe was calibrated with a Z-wire phantom by the manufacturer. Finally, the error due to rigid registration performed by the neuronavigation system has been estimated to be around 1.5mm in [START_REF] Lunn | Displacement Estimation with Co-Registered Ultrasound for Image Guided Neurosurgery: A Quantitative In Vivo Porcine[END_REF].

Pathology of the patients

In this study, hyperechogenic pathologies such as cavernoma (patient 1, see Fig. 5 and patient 2, see Fig. 6) and low-grade glioma (patient 3, see Fig. 7) were considered. In T1-w MR images, the central part of cavernoma is usually heterogeneous (hyper and hyposignal) and the outlying area appears in hyposignal. The low-grade glioma are more homogeneous and appears in hyposignal in T1-w MR images. In US images, numerous studies showed that all solid brain tumors, metastatic brain lesions, and cavernomas exhibited echogenicity [START_REF] Enzmann | Tumors of the central nervous system studied by computed tomography and ultrasound[END_REF][START_REF] Hammoud | Use of intraoperative ultrasound for localizing tumors and determining the extent of resection: a comparative study with magnetic resonance imaging[END_REF][START_REF] Le Roux | Low grade gliomas: comparison of intraoperative ultrasound characteristics with preoperative imaging studies[END_REF][START_REF] Mcgahan | Brain gliomas: sonographic characterization[END_REF][START_REF] Unsgaard | Ability of navigated 3d[END_REF]. For brain gliomas, the higher its grade (more malignant), the more echogenic it is in US and the less homogeneous it appears. In our study, the corresponding lesional tissues were considered both homogeneous and hyperechogenic in US images. As such, Ψ(X) was set to 1 for all segmentation of pathological tissue M 2 (see Eq. 5). Typical examples of intraoperative images and probability maps are presented in Fig. 2, Fig. 3 and Fig. 4.

Parameter Settings

The maximization of the joint probability (see Eq. 2) is performed within a multiresolution procedure using the simplex algorithm [START_REF] Nelder | A Simplex Method for Function Minimization[END_REF]. During the experiments, the parameters of the simplex algorithm were: tolerance = 0.1, stepsize = 1.5 and maximum number of iterations = 100. The coarsest resolution corresponded to the original volumes downsampled by a factor 3 and the finest resolution was that of the original volumes. The registration procedures take less than two minutes on Intel Pentium M at 2GHz.

As with most derivative-based operators, the M Lvv operator uses Gaussian kernel to compute the image derivatives. In [START_REF] Luc | Scale and the differential structure of images[END_REF], the authors showed that the convolution of the image with a derivative Gaussian kernel provides a wellposed approach of the differentiation problem. The standard deviation σ of the Gaussian kernel is called the image scale. This parameter has been shown as very stable for MR image sulci segmentation on numerous works [START_REF] Goualher | Statistical sulcal shape comparisons: Application to the detection of genetic encoding of the central sulcus shape[END_REF][START_REF] Goualher | Three-dimensional segmentation and representation of cortical sulci using active ribbons[END_REF][START_REF] Goualher | Automated extraction and variability analysis of sulcal neuroanatomy[END_REF], and thus, no tuning has been done for this parameter throughout our study. In our experiments, an image scale of 2 voxels has been used to compute the M Lvv values. This value is consistent with others works [START_REF] Goualher | Statistical sulcal shape comparisons: Application to the detection of genetic encoding of the central sulcus shape[END_REF][START_REF] Goualher | Three-dimensional segmentation and representation of cortical sulci using active ribbons[END_REF][START_REF] Goualher | Automated extraction and variability analysis of sulcal neuroanatomy[END_REF] conducted on brain cortical segmentation where the scale parameter was always kept in this range.

Evaluation Framework

In order to evaluate our method, a validation framework with different approaches is proposed.

• First, a visual assessment is proposed.

• Second, a manual validation by experts is presented. This validation is divided in two parts: a point-based estimation of the rigid registration by 3 experts for the 3 patients, and an evaluation of the residual error by all experts for 1 patient (post-registration error).

• Third, a study on convergence robustness was carried out.

The expert manual validation was difficult due to the time required. For each expert, 4 hours were required to perform the a priori estimation of the transformation for 3 patients.

Visual Assessment

The visual assessment remains a valuable indicator of the registration accuracy. In [START_REF] Holden | Voxel similarity measures for 3-D serial MR brain image registration[END_REF], the observer discernibility of registration errors has been estimated around 0.2 mm. A study on visual inspection for image registration assessment can be found in [START_REF] Fitzpatrick | Visual assessment of the accuracy of retrospective registration of MR and CT images of the brain[END_REF]. In our paper, we propose an overlay of US and MR images before and after registration to assess the registration accuracy.

Validation by experts

First, the experts manually evaluate the rigid transformation between the intraoperative US and the MR image resliced with the rigid transformation given by the neuronavigation system. This estimation is denoted as a priori estimation of the registration. From this a priori estimation, the initial error (i.e. after the registration performed by the neuronavigation system) and the Target Registration Error (TRE) can be computed. The a priori estimation of the registration is used to show that there are no statistical differences between the expert-based transformations and the transformation estimated with our method in terms of the TRE.

The experts estimate the residual error after rigid registration based on a given transformation (either by our method or the point-based expert registrations). This estimation is called a posteriori evaluation of the residual error and is designed to show that experts do not detect significant differences when they inspect the registered volumes with our method or with their own manually defined transformations.

A priori estimation of the registration

Point picking The a priori estimation of the registration is based on the location of ten points in the US image and the ten corresponding points in the MR image: each expert defines a set of point in the 3D reconstruction of the intraoperative ultrasound and its corresponding landmark in the resliced MR image. The resliced MR image is obtained with the rigid registration given by the neuronavigation system and has the same resolutions, dimensions, and field of view as the reconstructed US image. During the experiments, the experts used three orthogonal 2D views to define homologous points in the 3D volumes. For each volume, the visualization software was run independently, with the cursors in the two volumes unlinked. Each expert was allowed to choose their set of homologous points.

Initial error

The initial error is computed by using the mean Euclidean distance between the homologous points defined by the experts in both modalities. The three samples (one per expert) containing the ten error values (one per point) are compared by using a Kruskal-Wallis test.

Target Registration Error A leave-one-out procedure is used to compute the TRE of each point (i.e. Euclidean distance between homologous point after rigid transformation). First, one of the ten homologous point is removed from the set of points. The nine remaining homologous points are then used to compute a rigid transformation in the least squares sense. Finally, this rigid transformation is used to compute the TRE of the initially removed point. This procedure is repeated for all the ten points. The final TRE is the mean TRE over all the points. For each patient, the expert-based TRE and the TRE obtained with our method are compared by using a non parametric Kruskal-Wallis test.

A posteriori evaluation of the residual error

Point picking First, the patient images are registered using several transformations. These transformations are: i ) the three expert-based transformations ( T1 , T2 , T3 ), ii ) the rigid transformation obtained with our method ( T ), and iii ) the transformation computed using all the thirty points defined by the three experts Tall . The experts then define ten homologous points on the registered volumes. This procedure is performed for the five studied registrations on the patient 2 dataset. The positions of the points are fixed for all the experts.

Residual error As for initial error, the final error or residual error is simply obtained by computing the mean Euclidean distance between the homologous points defined by the experts in both modalities. The statistical comparison of the residual errors is performed on the five samples (one per transformation) containing ten errors values (one per point) with a Kruskal-Wallis test.

Robustness study

First, the US and resliced MR images of patient 2 are registered with the transformation Tall . Then, 100 rigid transformations are randomly generated with a translation along each axis uniformly distributed between 0 -5 mm and with a rotation around each axis uniformly distributed between 0 -5 degrees. Finally, each transformation is applied to the resliced MR image before performing registrations with the proposed method. The warping index ω [START_REF] Thévenaz | Optimization of mutual information for multiresolution image registration[END_REF] is used to compute the distance between the estimated transformation by the registration process T and the true transformation T :

ω = 1 |Ω| X∈Ω T -1 (X) -T (X) 2 (7)
where ||.|| 2 is the L 2 -norm. The success rate is estimated by considering a success as a registration with a warping index inferior to 3.5 mm. This threshold has been chosen close to the upper bound of the TRE estimated by the experts (see distribution for patient 2 in Fig. 9). Contrary to TRE estimated over selected points, the warping index is computed as the average error between the volumes over all the voxels.

3 Results

Visual assessment

The registration results are first displayed for visual assessment. The results obtained with our method are presented in Fig 5,6 and 7. For patient 1 (see Fig. 5), even if the lesion was not entirely included in the US volume, the proposed registration procedure converged efficiently. For patient 2 (see Fig. 6), acoustic shadows are present on the US image. The signal below the lesion tends to zero. The proposed approach overcomes these artifacts without specific detection of the shadows. For patient 3 (see Fig. 7), despite the large size of the low-grade glioma and the limited field of view, our approach performed well.

Validation by experts

A priori estimation of the registration

Table 1 presents the estimated initial error for the three patients by the three experts. The p-value of the Kruskal-Wallis test showed that there was no significant difference between the expert estimations. Table 1 also shows the interindividual variability for the same measure between the experts. The Fig. 8 summarizes the distribution of the error.

The estimated initial errors are significantly higher than values given in [5,28] (< 3mm) or by the manufacturer (< 1.5mm). It is important to note that the ultrasound images used in our experiments were acquired in clinical context during a neurosurgical operation. The real neurosurgery context is likely more difficult than phantom and animal studies. 1: Manual estimation of the initial error in mm (i.e. error of the registration given by the neuronavigation system).

A
Tab. 2 shows the TRE estimated by each expert, for each patient dataset. In all the cases, there were no statistically significant differences between the TRE obtained with expert-based estimations and the TRE obtained with our method. Figure 9 shows the result of the Kruskal-Wallis test. In all the cases, the experts and our method provided consistent results. Left : registration given by the neuronavigation system. Right : the result after correction with our registration approach. The low intensities of US images are in green and the high in red. For this case, even if the lesion was not entirely included in the US volume, the proposed registration procedure converged. The location of the lesion is indicated by ellipses on the fused image after registration.

x-y plane x-z plane y-z plane Figure 6: Patient 2. Left : registration given by the neuronavigation system. Right : results after correction with our registration approach. The low intensities of US images are in green and the high in red. In this case, the acoustic shadow artifact was present on the US image. The signal below the lesion was totally dark as shown on Fig. 3. The proposed approach allowed to overcome these artifacts without specific detection of the shadows.

x-y plane x-z plane y-z plane A priori estimation of the registration:

A posteriori evaluation of the residual error

Table 3 shows the expert-based estimation of the a posteriori residual error of the different registrations (manual-based T and automatic T ) proposed for patient 2. The Kruskal-Wallis test shows that the errors associated with the transformations ( T1 , T2 , T3 , Tall ) and T are not significantly different. Fig. 10 shows the statistical distribution of the residual error for each transformation compared. Finally, the experts failed to detect significant differences between the manual-based registrations and our automatic registration. The residual error estimated by experts is around 1 -1.5 mm for all the transformations.

Discusion and Conclusion

This paper presents a new framework for the 3D rigid registration of US and T1-w MR brain images. In order to address this changeling problem, we propose an innovative probabilistic objective function that maximizes the joint probability of the i) a priori most probable locations of hyperechogenic structure in the preoperative MR image, and ii ) the highest intensities of the intraoperative US images. We show that the proposed method enables a robust registration of MR and US images in a computational time compatible with clinical use. All our experiments were carried out on real intraoperative data. The expert-based quantitative study shows that our method produces no statistically different registration compared to the a priori estimation of the registration by the experts. Moreover, the a posteriori estimation of the residual registration error shows that the experts failed to detect differences between manual registration and our automatic registration.

During our experiments, manual segmentation has been used to built the probability map. This segmentation is always available, since the neurosurgeon performed it before the surgery. In this paper, the used segmentations were the segmentations dedicated to the neurosurgery. However, the segmentation of the lesion could be automated [START_REF] Kaus | Automated segmentation of MR images of brain tumors[END_REF][START_REF] Prastawa | Automatic brain tumor segmentation by subject specific modification of atlas priors[END_REF] and the different parts of pathologies (lesion, coagulated blood, cyst, necrotic tissue, etc) could be defined. Through this, the simple model of homogeneous hyperechogenic lesion used in our experiment could be improved by using different hyperechogenic levels to the different pathological tissues. To evaluate the robustness of our method to heterogeneous lesion, more datasets are needed although this situation was present in case of Patient 1.

The proposed method is related to the segmentation accuracy of the tumor in preoperative MR images. Although the segmentation of the MR image is not considered as difficult, this step may introduce some errors. Our experiments showed that the proposed method produced consistent results with manual seg-mentation used in clinical routine.

The presented clinical datasets showed that our method is robust to some discrepancies between the features present in both US and MRI probability maps. Based on the correlation of maps where only regions considered as relevant are used to drive the registration procedure, our methods is able to deal with partially missing information resulting from a limited field of view or acoustic shadows. In case of Patient 1 (see Fig 2 and Fig. 5), only a subpart of the lesion was visible in the reconstructed US image. In case of Patient 2 (see Fig. 3 and Fig. 6) the acoustic shadow below the lesion reduced information around ventricle in US. Moreover, the information derived from sulci was much more present in MR maps than in US map. Finally, in case of Patient 3 (see Fig. 4 and Fig. 7) the limited field of view and the large size of glioma reduced the importance of sulcal information. However, experiments using only the segmentation of the lesion or only sulcal information derived from M lvv operator failed to provide satisfactory registration. This seems to indicate that a certain amount of homologous features has to be present in both probability maps to enable the method working.

Finally, in our opinion, the proposed approach relies on a similar and complementary idea to the vessel-based method proposed by Reinertsen's et al. [START_REF] Reinertsen | Validation of vessel-based registration for correction of brain shift[END_REF][START_REF] Reinertsen | Clinical validation of vessel-based registration for correction of brain-shift[END_REF]. Indeed, in both cases, an implicit segmentation of salient features in US images (hyperechogenic structures in B-mode or vessels in Doppler) is matched with corresponding structures detected in MR images. Only the selected salient features differ between the methods. In [START_REF] Reinertsen | Validation of vessel-based registration for correction of brain shift[END_REF][START_REF] Reinertsen | Clinical validation of vessel-based registration for correction of brain-shift[END_REF], the method utilizes vessels extracted from Doppler US images and their segmentation from MR images. Therefore, both methods have the advantage of not requiring segmentation of the US image and also being robust to US artefacts. However, the extraction of the vessel centerlines from MR images is a challenging problem and requires extensive processing.

Our method is dedicated to brain US imaging since M Lvv operator is relevant for sulci and cerebral falx detection. As such, the application of the proposed framework to another body part requires adaptation of the hyperechogenic structure detection. Moreover, if T2-w MR image or another sequence is used as preoperative MR image, the selected values of the M Lvv need to be adapted. Since the final aim of this US/MR registration method is to compensate for the brainshift, further works will investigate extension of our probabilistic objective function to non-rigid deformations. ultrasound to delineate gliomas and metastases-comparison of image interpretations with histopathology. Acta neurochirurgica, 147 [START_REF] Luc | Scale and the differential structure of images[END_REF], December 2005.

Figure 1 :

 1 Figure 1: Illustration of the performed workflow to achieve the registration. The skull stripping, the denoising, the M Lvv computation and the segmentation of lesion are performed before the neurosurgical procedure. The 3D reconstruction of intraoperative volume, the reslicing of the MR map and the estimation of the transformation are then performed during the neurosurgical procedure.

Figure 2 :

 2 Figure 2: Patient 1.Left: the denoised MR image corresponding to US image (Middle). Right: the probability map based on M Lvv operator and extracted from denoising MR image. The preoperative MR image is resliced with the registration matrix provided by the neuronavigation system. The matching is not perfect due to the initial error, estimated around 4.2 mm by the experts for this patient. As visible on the US image, the lesion presents a hypoechogenic central area (indicated on the US image by a narrow).

Figure 3 :

 3 Figure 3: Patient 2. Left: the denoised MR image corresponding to US image (Middle). Right: the probability map based on M Lvv operator and extracted from denoising MR image. The preoperative MR image is resliced with the registration matrix provided by the neuronavigation system. The matching is not perfect due to the initial error, estimated around 8.5 mm by the experts for this patient. The red box shows a sulcal area where the matching between US image and the M Lvv-based probability map is visually high. A large acoustic shadow is visible in US around the ventricle area (indicated in orange on the US image).

Figure 4 :Figure 5 :

 45 Figure 4: Patient 3. Left: the denoised MR image corresponding to US image (Middle). Right: the probability map based on M Lvv operator and extracted from denoising MR image. The initial error is estimated around 5 mm by the experts for this patient. In the red box, the M Lvv operator efficiently detects the sulci also visible in the US image. This case presents a limited field of view while the lesion is large. Moreover, an acoustic shadow is visible on the right lower part of the US image ( (indicated on the US image by a narrow).)

Figure 7 :Figure 8 :

 78 Figure 7: Patient 3. Left : registration given by the neuronavigation system. Right : the result after correction with our registration approach.The low intensities of US images are in green and the high in red.

Figure 9 :

 9 Figure 9: Target Registration Error (in mm) estimated by all experts for the each patient dataset.

Figure 11 :

 11 Figure 11: Distribution of the warping index in mm obtained with our method.

  priori estimation of the registration: Initial error in mm

	mean (std)	Expert 1	Expert 2	Expert 3	p-value
	Patient 1	5.52 (1.15) 4.31 (1.55) 5.00 (1.50)	0.30
	Patient 2	8.64 (0.89) 8.31 (1.24) 8.76 (1.04)	0.69
	Patient 3	3.56 (1.09) 4.61 (1.39) 4.38 (1.13)	0.13
	Table

Table 2 :

 2 Target Registration Error in mm. The p-values correspond to Kruskal-Wallis test perform between the TRE obtained by experts and the TRE obtained with our method.

	Target Registration Error in mm

https://www.irisa.fr/visages/benchmarks/

Acknowledgment

We are grateful for Sean Jy-Shyang Chen for proofreading the manuscript.

Robustness study

Tab. [START_REF] Cachia | A generic framework for parcellation of the cortical surface into gyri using geodesic Voronoi diagrams[END_REF] shows the robustness and the warping index results obtained during the experiment. The proposed method obtained 92% of success rate with a mean warping index of 2.38 mm. This value is relative to the TRE of the used gold standard. Thus, it gives information about the distance between the transformation from all the experts ( Tall ) and the final transformation provided by our method. This value is close to the TRE estimated for patient 2 in Tab.

2. The figure 11 shows the distribution of the warping index.

Success rate in % ω in mm (mean (std))