B Dufay

S Saez

C Dolabdjian

A Yelon

D Ménard

Impact of electronic conditioning on the noise performance of a two-port network Giant MagnetoImpedance magnetometer

Keywords: Magnetometer, GMI, Noise

The performance of GMI-based magnetometers is currently limited by the noise due to the electronic conditioning circuitry. We propose a simple model of this noise for a GMI sensor using a synchronous detection scheme. The GMI sensing element consists of a thin pick-up coil wound around a Co-rich amorphous micro-wire. It is fully described by a two port network model and associated impedance matrix. Noise and sensitivity behavior are studied for the four measuring configurations, corresponding to four terms of the impedance matrix. The model yields a good description of experimental data from noise measurements. The magnetic noise spectral density is dominated either by the excitation or detection stages, depending upon whether the excitation currents are high or low. The non-trivial noise behavior exhibited by each configuration leads to better understanding of the noise limitations of GMI magnetometers. The configuration in which the signal at the coil terminals is measured (often named off-diagonal) is the most efficient in decreasing the equivalent output magnetic noise spectral density.

I. INTRODUCTION

The Giant Magneto-Impedance (GMI) effect reflects the strong variation of the impedance of a ferromagnetic conductor, driven by a high frequency current, as a function of an external applied magnetic field parallel to the wire. This strong effect, which is very promising for the development of highly sensitive magnetometers [START_REF] Knobel | Giant magnetoimpedance[END_REF], [START_REF] Bushida | Sensitive and quick response micro magnetic sensor using amorphous wire MI element colpitts oscillator[END_REF], [START_REF] Hauser | Giant magnetoimpedance sensors[END_REF], [START_REF] Mohri | Amorphous wire and cmos ic-based sensitive micromagnetic sensors utilizing magnetoimpedance (mi) and stressimpedance (si) effects[END_REF], [START_REF] Yabukami | A thin film magnetic field sensor of sub-pt resolution and magnetocardiogram (mcg) measurement at room temperature[END_REF], has thus attracted considerable attention over the past few years. A very important characteristic of magnetometers is their equivalent magnetic noise level, which ultimately limits the sensing performance. We previously reported that the noise level of current GMI-based magnetometers, in a typical measurement setup, is dominated by the noise of the electronic conditioning, rather than by the intrinsic noise of the GMI sensing element [START_REF] Melo | Optimization of the magnetic noise and sensitivity of giant magnetoimpedance sensors[END_REF], [START_REF] Ding | Low frequency giant magnetoimpedance magnetometer noise versus electronic conditioning[END_REF], [START_REF] Ding | Investigation of giant magnetoimpedance magnetic noise comparison[END_REF]. As a consequence, it is crucial to increase the intrinsic sensitivity of the sensor ∂Z ∂B (i.e. the differential variation of the impedance expressed in Ω/T) until the sensor noise becomes higher than the electronic conditioning noise, in order to further reduce the total equivalent magnetic noise. With that goal in mind, we investigated the noise of a complete GMI device (coupled with thin pick-up coil, or not). Our approach is based on previous work, in which we showed that the GMI-coil set-up, generally labeled off-diagonal GMI [START_REF] Sandacci | Off-diagonal impedance in amorphous wires and its application to linear magnetic sensors[END_REF], [START_REF] Kraus | Off-diagonal magnetoimpedance in stress-annealed amorphous ribbons[END_REF], orthogonal flux-gate in the fundamental mode [START_REF] Zhao | Comparative study of the sensing performance of orthogonal fluxgate sensors with different amorphous sensing elements[END_REF], [START_REF] Sasada | Orthogonal fluxgate mechanism operated with dc biased excitation[END_REF], or inverse Wiedemann effect [START_REF] Kraus | Magnetic field sensor based on asymmetric inverse Wiedemann effect[END_REF], is a promising approach to increasing the intrinsic sensitivity [START_REF] Dufay | Improved GMI sensors using strongly-coupled thin pick-up coil[END_REF]. Here, more attention is dedicated to the noise aspect of the detection configuration, taking into account the electronic conditioning limitation, evaluating the impact of using or not using the coupled pick-up coil. The focus of this article is thus to present the expected noise behavior of a GMI device, including the effect of electronic conditioning. Actual optimization of a complete device, based on the result of this study, will be presented elsewhere.

The paper is organized as follows. Section II describes the electronic setup. The associated noise model is developped in Sec. III. Preliminary results and discussion are presented in section IV. Methods of noise optimisation are discussed in section V, which is followed by a general conclusion.

II. EXPERIMENTAL METHOD

Here we present the electronic setup used in our study. We first describe the sensor itself, which consists of a GMI micro-wire coupled to a thin pick-up coil. We then present the electronic conditioning circuit which sets the sensor in appropriate operating conditions.

The sensing element consists of a very thin bi-layer pick-up coil, with approximately 500 turns per layer, directly wound directly around a 100 µm diameter CoFeSiB amorphous ferromagnetic wire. The length of the complete device is approximately 2.5 cm. A schematic of the device along with its equivalent two-port network model are shown in Fig. 1.

The scattering parameters of the device were measured using a classical test setup based on a vector network anal-yser [START_REF] Dufay | Improved GMI sensors using strongly-coupled thin pick-up coil[END_REF], driven with sufficiently low excitation power as to insure that the GMI element operated in the linear regime. A description in terms of a field-dependent impedance matrix [Z(B)], where B (B = µ 0 H) is the external magnetic induction, links the voltages to the currents (see Fig. 1 for definitions), that is

v 1 v 2 = [Z(B)] i 1 i 2 = Z 11 (B) Z 12 (B) Z 21 (B) Z 22 (B) i 1 i 2 . ( 1
)
Typical results of measurements of the components of this impedance matrix are illustrated in Fig. 2. Several dc bias currents are superposed on the excitation current, of frequency, f 0 = 300 kHz. This relatively low driving frequency has been chosen in order to minimize coupling (crosstalk) due to wave propagation, despite the improvement in performance expected at higher frequencies. Our focus here is on comprehension of the noise behavior rather than on the optimization of the response. Let us consider that the circuit is operated in a fieldlocked loop configuration. The external magnetic field B is written as

B = B 0 + b(t) (2)
where B 0 is the static working point and b(t) is the small magnetic field signal [START_REF] Ding | Equivalent magnetic noise limit of low-cost GMI magnetometer[END_REF]. A first order Taylor expansion allows us to write each term of the impedance matrix as

Z ij (B) = Z ij0 + ∂Z ij (B) ∂B B=B0 .b(t) (3)
where

Z ij0 = Z ij (B 0 ) and ∂Zij (B) ∂B B=B0
is the differential impedance variation as a function of the applied magnetic field, expressed in units of Ω/T. As suggested by the curves of Fig. 2, a dc bias current of +6 mA is chosen for the following discussion, in order to enhance the differential impedance variation [START_REF] Melo | Optimization of the magnetic noise and sensitivity of giant magnetoimpedance sensors[END_REF].

In addition to the sensing element, the electronic conditioning circuit consists of an excitation stage, providing an ac current to the two-port network, and a detection stage, which collects the output signal. The schematic of this circuit is shown in Fig. 3. The excitation stage consists of a 300 kHz sine wave voltage generator (HP 33120A), having an amplitude, e g1 , of 1 V p . This constant voltage source acts as a constant current source, due to the resistor R 1 . This equivalent sine wave current, which flows through the twoport network, induces a sine wave voltage at output ports. This voltage is modulated in amplitude, with a carrier frequency f 0 , where the modulation depth is directly related to the terms ∂Zij ∂B B=B0

.b(t) of the impedance matrix.

The detection stage consists of a preamplifier (EG&G Model 5113) which picks up the voltage at the two-port network terminals. Its use is motivated by the need for high input impedance compared to each Z ij0 . A large gain, G ≈ 1000, is chosen, allowing us to neglect any noise sources after the preamplifier. While the proposed model is linear, the experimental behavior for very small R 1 may not be. To minimize this problem, a second order low-pass filter with a cut-off frequency of f 0 , follows the preamplifier, in order to suppress any small distortion appearing at high excitation current amplitudes due to non-linear effects in the GMI wire. In order to avoid saturation of the preamplifier and because of this large gain G, we need to cancel out the constant average amplitude of the carrier, keeping only the modulated part of the amplitude. This explains the second generator voltage connected to the inverting input of the preamplifier, synchronized to the first, and at the same frequency, with a voltage amplitude, e g2 .

Following the preamplifier, we use a Stanford Research System (SR 844) Lock-In Amplifier (LIA), locked to the carrier frequency, f 0 . It provides an output signal, directly proportional to the modulation amplitude of the input signal, with gain, G LIA . From this description of the system, we have developed a noise model, presented in the next section.

III. MODELING OF NOISE SOURCES AND EVALUATION OF THE SENSITIVITY

The noise model developed here allows one to predict noise performance of the setup described in the previous section. First, we quantify the contribution of each noise source, relative to the white noise level at the output of the setup after demodulation. We focus on the voltage spectral density (white noise) in a range around ±1 to 10 kHz from the carrier, f 0 , which will be shifted to low frequency after demodulation. Then we predict the output sensitivity, T r , in V/T. This links the output signal, V out (t), to the applied external magnetic field, b(t).

In all of the following, we identify the different excitation/detection configurations by a subscript X (X = A, B, C, or D), as defined in Fig. 3 and referring to the impedance matrix term given in Fig. 1. In brief, A means excitation through the amorphous wire with signal detection across its terminals, B means the same excitation, but with detection across the coil terminals. C means the opposite of B. D means excitation through the coil combined with detection across its terminals.

Assuming that the intrinsic noise level of the sensing element is lower than that of the conditioning electronic setup [START_REF] Ding | Intrinsic giant magnetoimpedance noise reduction by dc bias[END_REF], we identify three kinds of noise sources. The first is noise induced by the voltage generators, e g . Signal instability Figure 2. Real and imaginary parts and modulus of Z ij terms of the impedance matrix as a function of the applied magnetic field for three dc bias currents flowing through the amorphous wire. Measurements were performed for an excitation frequency, f 0 , of 300 kHz. On the modulus curves, we show estimated differential variations of impedance expressed in Ω/T. These are used in section III to estimate the sensitivity. of sinusoidal sources is generally characterized by the singlesideband noise spectral density (NSD) expressed in decibel below the carrier per hertz (dB c /Hz) [START_REF] Lee | Oscillator phase noise: a tutorial[END_REF], in direct relation to the output amplitude of the source. In our case, we consider a dynamic range of 148 dBc/Hz [START_REF]Function / Arbitrary waveform generator HP 33120A User's Guide[END_REF], yielding the voltage power NSD of the generator, expressed as

e 2 ng1,2 = eg 1,2 √ 2 2 10 148/10 = e g1,2 10 151/20 2 [V 2 rms /Hz] (4)
where e g1,2 is the amplitude of generator 1 or 2 in V peak . The amplitude, e g2 , of the second generator is related to e g1 since amplitude of the signals at both inverting and non-inverting inputs of the preamplifier need to be approximately equal.

Consequently, the noise level of the second generator, e ng2 , may be expressed as functions of e ng1 and of circuit elements.

The second noise source to be taken into account is the noise level of the preamplifier, which may be summarized by its (e n , i n ) model, considering an input voltage white power

e 2 nR = 4k B T R [V 2 rms /Hz] (5)
where k B ≈ 1.38 × 10 -23 JK -1 is the Boltzmann constant and T ≈ 300 K is the operating temperature. Assuming that the input impedance of the preamplifier is sufficiently high that we may neglect input current and that the various noise sources are uncorrelated, the white NSD close to the carrier frequency (∆f ≈ ±1 to 10 kHz) at the LIA input can be obtained from Eqs. (6a) to (6d), for each configuration, X:

In these expressions, noise due to internal generator resistors (50 Ω) is neglected compared to the generator noise level itself.

e 2 nLIAA (f 0 ± ∆f ) ≈ G 2 |Z11 0 (f0)| 2 (R1+50+|Z11 0 (f0)|) 2 2e 2 ng1 + e 2 nR1 + (R 1 + 50) 2 i 2 npreamp +e 2 npreamp + e 2 nR=1000 + 1000 2 i 2 npreamp , (6a)
e 2 nLIAB (f 0 ± ∆f ) ≈ G 2 |Z21 0 (f0)| 2 (R1+50+|Z11 0 (f0)|) 2 2e 2 ng1 + e 2 nR1 + e 2 nR=1000 + |Z 220 (f 0 )| - |Z21 0 (f0)||Z12 0 (f0)| R1+50+|Z11 0 (f0)| 2 i 2 npreamp + e 2 npreamp + 1000 2 i 2 npreamp , (6b)
e 2 nLIAC (f 0 ± ∆f ) ≈ G 2 |Z12 0 (f0)| 2 (R1+50+|Z22 0 (f0)|) 2 2e 2 ng1 + e 2 nR1 + e 2 nR=1000 + |Z 110 (f 0 )| - |Z21 0 (f0)||Z12 0 (f0)| R1+50+|Z22 0 (f0)| 2 i 2 npreamp + e 2 npreamp + 1000 2 i 2 npreamp , (6c)
e 2 nLIAD (f 0 ± ∆f ) ≈ G 2 |Z22 0 (f0)| 2 (R1+50+|Z22 0 (f0)|) 2 2e 2 ng1 + e 2 nR1 + (R 1 + 50) 2 i 2 npreamp +e 2 npreamp + e 2 nR=1000 + 1000 2 i 2 npreamp . (6d)
Considering an AM signal at the LIA input of the form [1 + A cos(ωt)] cos(ω 0 t), where ω 0 and ω are the carrier and modulating pulsation respectively, the filtered demodulated signal is G LIA A cos(ωt), where G LIA is related to the gain of the demodulation stage. This output signal is due to the sum of both upper and lower sidebands, each sideband being correlated with the same amplitude A 2 . However, for an input white noise around the carrier, upper and lower sidebands are not correlated. Consequently, the output NSD is increased by a factor G LIA √ 2, due to the quadratic sum of both the noise of the two sidebands. This effectively results in a decrease by a factor √ 2 of the signal to noise ratio 1 . Consequently, we express the equivalent output white noise power spectral density given at the sensor output by

e 2 nX (∆f ) = √ 2G LIA G LIA G 2 e 2 nLIAX (f 0 ± ∆f ) = √ 2 G 2 e 2 nLIAX (f 0 ± ∆f ). (7)
The next step is to evaluate the output sensitivity at the sensor output, represented by

T rX = ∂Vout X (t) ∂B 1
GLIAG , in V/T. This is easily evaluated using ∂Zij ∂B and Z ij0 of Eq. 3, and from circuit elements. Equations 8a to 8d give the sensitivity of each configuration (A, B, C, or D).

Finally, one can evaluate the equivalent magnetic noise spectral density of the setup, b nX , in pT/ √ Hz. It is defined as the ratio of the electronic noise spectral density (in V/ √ Hz) to the sensitivity (in V/T), b nX = en X Tr X .

IV. PRELIMINARY RESULTS AND DISCUSSION

With the help of the simplified model, developed in the previous section, we can now analyze the noise behavior for the various excitation conditions, which is determined by the injection resistance, R 1 , and the generator voltage amplitude, e g1 . In this section, we consider a static working point, B 0 , which sets static elements of the impedance matrix, Z ij0 , to half maximum where differential variations of impedance,

∂Zij (B) ∂B B=B0
, are maximum. This is illustrated by red arrows on modulus curves of Fig. 2. Fig. 4 shows all noise contributions to the total NSD before demodulation, e nLIAX , for all four configurations as a function of R 1 and for an excitation e g1 = 1 V peak . We may easily identify the dominant noise sources for each excitation condition. As expected, the NSD is dominated by generator noise for high excitation currents (R 1 < 1 kΩ for e nA ) and by detection stage noise for lower excitation currents (combination of both preamplifier voltage noise, e npreamp , and 1 kΩ resistor noise). Equivalent input voltage NSDs after demodulation, e nX , of each configuration, are summarized in Fig. 5.

Figure 6 shows the sensitivity versus excitation level for all four configurations. Note that the very high intrinsic sensitivity ∂Z22 ∂B , shown in Fig. 2, does not necessarily lead to the highest output sensitivity for low values of R 1 because of the high value of Z 220 . This, among other things, illustrates the necessity of fully understanding the entire system in order to choose the optimal conditioning for each configuration.

Finally, Fig. 7 shows the calculated equivalent magnetic noise levels in pT/ √ Hz as a function of R 1 . The figure shows non-trivial behavior of the noise with respect to electronic conditioning. It appears that the B configuration (GMI → Coil) is the most promising for the set of parameters considered here. The optimal operating point (excitation) is located in the lowest flat part of the curve, where the dominating noise is due to the excitation stage. It is important to point out that, for an equivalent white NSD dominated by excitation noise, it is useless to further increase the excitation amplitude (in the linear regime) to produce higher response, because the noise level will increase as the sensitivity does. For the lowest values of R 1 , the current flowing through the sample , at ∆f = ±10 kHz from the carrier, f 0 , expressed at the sensor output (a) for each configuration, X.

T rA (f 0 ) ≈ ∂ |Z 11 (f 0 )| ∂B e g1 (R 1 + 50) (|Z 110 (f 0 )| + R 1 + 50) 2 , (8a)
T rB (f 0 ) ≈ e g1 |Z 110 (f 0 )| + R 1 + 50 ∂ |Z 21 (f 0 )| ∂B - |Z 210 (f 0 )| |Z 110 (f 0 )| + R 1 + 50 ∂ |Z 11 (f 0 )| ∂B (8b) ≈ e g1 |Z 110 (f 0 )| + R 1 + 50 ∂ |Z 21 (f 0 )| ∂B , T rC (f 0 ) ≈ e g1 |Z 220 (f 0 )| + R 1 + 50 ∂ |Z 12 (f 0 )| ∂B - |Z 120 (f 0 )| |Z 220 (f 0 )| + R 1 + 50 ∂ |Z 22 (f 0 )| ∂B (8c) ≈ e g1 |Z 220 (f 0 )| + R 1 + 50 ∂ |Z 12 (f 0 )| ∂B , T rD (f 0 ) ≈ ∂ |Z 22 (f 0 )| ∂B e g1 (R 1 + 50) (|Z 220 (f 0 )| + R 1 + 50) 2 . (8d)

Curves (b)-(f) indicate the various noise contributions to

e nLIA X G
as determined from Eqs. (6a) to (6d). These are the noise arising from the 1 kΩ resistor and network elements, (b), from the preamplifier current noise and network elements, (c), from the preamplifier voltage noise and network elements, (d), from generators and network elements, (e), and from the resistor R 1 and network elements, (f). These noise levels are shown versus the injection resistance R 1 for an excitation amplitude, eg 1 , of 1 V peak . Figure 5. Theoretical (line) and measured (symbols) equivalent white noise spectral density, expressed at the sensor output, in V/ √ Hz versus resistor value, R 1 , at an excitation amplitude, eg 1 , of 1 V peak . Each curve corresponds to one of the four available configuration cases (A, B, C, or D). Filled symbols are measured before LIA and considering an optimal demodulation as described in the present model whereas empty symbols are measured after SR844 LIA. In this later case, the noise spectral density could be limited by the dynamic of the used LIA. In all the following, we will consider experimental NSD measured as filled symbols. Theoretical (line) and measured (symbols) sensitivity after demodulation, expressed at the sensor output, in V/T versus resistor value, R 1 , at an excitation amplitude, eg 1 , of 1 V peak . Each curve corresponds to one of the four available configuration (A, B, C, or D).

becomes too large and leads to non-linear GMI behavior with some distortion. Nevertheless, previous work on non-linear GMI [START_REF] Seddaoui | Nonlinear electromagnetic response of ferromagnetic metals: Magnetoimpedance in microwires[END_REF] suggests that the sensitivity will be maximum when sensing element is operated in the low non-linear range, just before the sensitivity starts to decrease. As a consequence, the best signal to noise ratio may be obtained at higher excitation current. However, this is not considered in the present model.

The model proposed here has been confronted with a series of noise measurements carried out on CoFeSiB microwire in the four configurations A, B, C, and D. The experimental details are given in Section II. The results are indicated by the symbols for each curve of Figs. 5 to 7. The good agreement between measured data and the model suggests that the overall approach to noise modeling is sound and that the model includes the essential mechanisms limiting the performance of Theoretical (line) and measured (symbols) equivalent white magnetic noise spectral density at the sensor output, in T/ √ Hz versus resistor value, R 1 , at an excitation amplitude, eg 1 , of 1 V peak . Each curve corresponds to one of the four available configurations (A, B, C, or D).

our GMI sensor. It also confirms, a posteriori, our assumption that the noise performance of the proposed design is still limited by electronic conditioning.

V. NOISE IMPROVEMENT In the previous section, we considered a static working point, B 0 , at which the Z ij0 are equal to their half maximum values. This choice is appropriate for validation of the model and discussion of the four configurations. Nevertheless, there is a small field range around this working point, B 0 , in which the impedance variation is essentially linear. In other words, we may assume that the differential impedance variation, ∂Zij (B) ∂B B=B0

, remains constant in this range whereas Z ij0 varies with B 0 . From Eqs. 6a to 6d, we note that the values of Z ij0 affect the contribution of generators noise to the equivalent output NSD. These static values likewise affect the sensitivity as shown in Eqs. 8a to 8d.

Using the noise model developed here, we can evaluate the equivalent magnetic NSD at each limit of the linear range, defining a range of achievable values surrounding the curves of Fig 7. This is shown in Fig. 8.

Figure 8 clearly shows that improved noise performance may be obtained, in particular for the B configuration (GMI → Coil) for which a decrease in Z 210 leads to a decrease in generator noise, the dominant noise source at high excitation current. This is confirmed by experimental measurements, as shown in Fig. 8, which compares results for the half maximum static working point (squares) with a static working point with a lower value of Z 210 (circles). As expected, there is a significant decrease of the equivalent magnetic NSD for the lowest values of excitation resistor R 1 . The best measured value, around 1 pT/ √ Hz, is not as good as expected. Indeed, for highest excitation current we cannot reach a value of Z 210 as low as expected without decreasing the sensitivity. VI. CONCLUSION Following our earlier finding, that the GMI-coil configuration (Case B above) constitutes a promising approach to increasing the sensitivity of GMI effect sensors [START_REF] Dufay | Improved GMI sensors using strongly-coupled thin pick-up coil[END_REF], we have demonstrated here that the noise performance may also be improved with this configuration. Indeed, our model allows us to understand the noise behavior with respect to four electronic conditioning configurations.

The good agreement between the model and experimental data validates our general approach. It thus allows us to extrapolate the results to other conditioning approaches and to predict the expected performance. For the present setup, the model leads to the determination of the most appropriate conditioning, providing the best possible performance for a given sensor configuration. It also highlights the dominant noise sources and limitations of GMI magnetometers. In the current state of the art, the noise performance of GMI magnetometers appears to be limited by that of the exciting generators. White noise levels as low as 1 pT √ Hz have been obtained in this nonoptimal conditioning configuration, leaving room for further improvement.

Figure 1 .

 1 Figure 1. Schematic of the device and associated electrical two port network model. It shows the different terms related to the impedance matrix given in (1).

Figure 3 .

 3 Figure3. Schematic of the electronic conditioning circuit. The excitation stage uses a voltage generator, eg 1 , and a resistor, R 1 . The detection stage consists of a preamplifier and associated carrier compensation circuitry, followed by a lock-in amplifier locked to the excitation frequency, f 0 .

Figure 4 .

 4 Figure 4. White noise spectral density, e nLIA X G

Figure 6 .

 6 Figure 6.Theoretical (line) and measured (symbols) sensitivity after demodulation, expressed at the sensor output, in V/T versus resistor value, R 1 , at an excitation amplitude, eg 1 , of 1 V peak . Each curve corresponds to one of the four available configuration (A, B, C, or D).

Figure 7 .

 7 Figure 7.Theoretical (line) and measured (symbols) equivalent white magnetic noise spectral density at the sensor output, in T/ √ Hz versus resistor value, R 1 , at an excitation amplitude, eg 1 , of 1 V peak . Each curve corresponds to one of the four available configurations (A, B, C, or D).

Figure 8 .

 8 Figure 8. Theoretical (line and area) and measured (symbols) equivalent white magnetic noise spectral density expressed at the sensor output, in T/ √ Hz versus resistor value, R 1 , at an excitation amplitude, eg 1 , of 1 V peak . Each curve corresponds to one of the four available configurations (A, B, C, or D). Colored areas illustrate achievable values for the static working point range in which the differential impedance variation is constant. Square and circle symbols indicate measured values in configuration B, for two different static working points, B 0 . The squares correspond to the working point of Fig. 7, and the circles to a lower value of Z 21 0 .

  The final noise source is related to the Johnson noise of each resistor, R, of the setup, expressed as

	NSD, e 2 npreamp , of 3.2 nV rms / NSD, i 2 npreamp , lower than 20 fA rms / √ Hz 2 √ and an input current Hz

2

[START_REF]Low Noise Voltage Preamplifier EG&G model 5113 data sheet[END_REF]

.

In case of demodulation with high carrier amplitude, signal to noise ratio may be reduced by the dynamics of the LIA.

ACKNOWLEDGMENT

The amorphous wire, based on CoFeSiB alloy, was kindly provided by H. Chiriac from National Institute of Research & Development for Technical Physics, Romania.


Facundo Mémoli
email: memoli@math.stanford.edu

Guillermo Sapiro

Paul Thompson
email: thompson@loni.ucla.edu

Geometric Surface and Brain Warping via Geodesic Minimizing Lipschitz Extensions

Based on the notion of Minimizing Lipschitz Extensions and its connection with the infinity Laplacian, a theoretical and computational framework for geometric nonrigid surface warping, and in particular the nonlinear registration of brain imaging data, is presented in this paper. The basic concept is to compute a map between surfaces that minimizes a distortion measure based on geodesic distances while respecting the provided boundary conditions. In particular, the global Lipschitz constant of the map is minimized. This framework allows generic boundary conditions to be applied and direct surfaceto-surface warping. It avoids the need for intermediate maps that flatten the surface onto the plane or sphere, as is commonly done in the literature on surface-based non-rigid registration. The presentation of the framework is complemented with examples on synthetic geometric phantoms and cortical surfaces extracted from human brain MRI scans.

Introduction

Nonrigid surface warping is one of the most fundamental problems in geometric surface processing. This is particularly relevant for problems such as shape comparison, motion and deformation analysis, and shape morphing and interpolation. A particular important example is when the surface represents a brain, for example, the boundary between white matter and gray matter or between gray matter and CSF. Brain warping, which is a form of brain image registration and geometric pattern matching that will be used as test case in this work, is one of the most fundamental and thereby most studied problems in computational brain imaging [START_REF] Toga | Brain Warping[END_REF]. Brain images are commonly warped, using 3-dimensional deformation fields, onto a common neuroanatomic template prior to cross-subject comparison and integration of functional and anatomical data. Images of the same subject may be warped into correspondence over time, to help analyze shape changes during development or degenerative diseases. Almost all the active research groups in this area have developed and/or have their favorite brain warping technique. 4 A few representative works can be found at [START_REF] Christensen | A deformable neuroanatomy textbook based on viscous fluid mechanics[END_REF][START_REF] Chung | Deformation-based surface morphometry applied to gray matter deformation[END_REF][START_REF] Collins | An automated 3D non-linear image deformation procedure for determination of gross morphometric variability in the human brain[END_REF][START_REF] Fischl | High-resolution inter-subject averaging and a coordinate system for the cortical surface[END_REF][START_REF] Freeborough | Modeling brain deformations in Alzheimer's disease by fluid registration of serial 3D MR images[END_REF][START_REF] Gu | Genus zero surface conformal mapping and its application to brain surface mapping[END_REF][START_REF] Miller | On the metrics and Euler-Lagrange equations of computational anatomy[END_REF][START_REF] Thompson | Growth patterns in the developing brain detected by using continuum-mechanical tensor maps[END_REF][START_REF] Thompson | Mathematical/computational challenges in creating population-based brain atlases[END_REF][START_REF] Toga | Brain Warping[END_REF][START_REF] Van Essen | Comparisons between human and macaque using shapebased deformation algorithms applied to cortical flat maps[END_REF][START_REF] Van Essen | Functional and structural mapping of human cerebral cortex: Solutions are in the surfaces[END_REF], this list being far from complete. In spite of this, the problem is still open and widely studied, since there is not a "ground truth" method to obtain a map between brains. The same is true for generic nonrigid surface matching. The criteria for matching different features (e.g., geometry or intensity) may also depend on the application, which range from recovering intraoperative brain change to mapping brain growth, or reducing cross-subject anatomical differences in group functional MRI studies.

The way the brain warping problem is addressed is critical for studies of brain diseases that are based on population comparisons. Examples of this application can be found at [START_REF] Freeborough | Modeling brain deformations in Alzheimer's disease by fluid registration of serial 3D MR images[END_REF][START_REF] Thompson | Growth patterns in the developing brain detected by using continuum-mechanical tensor maps[END_REF], although these are a very non-exhaustive account of the rich literature on the subject. The interested reader may also check [42] for numerous applications of brain warping and population studies. As detailed in [START_REF] Toga | Brain Warping[END_REF], shape and brain warping approaches can be divided into two classes, those based on volume-to-volume matching and those based on surface-to-surface matching. Our work belongs to the latter of the categories. Surface matching has recently received increasing attention as most functional brain imaging studies focus on the cortex, which varies widely in geometry across subjects. The power of these studies depends on the degree to which the functional anatomy of the cortex can be aligned across subjects, so improved cortical surface registration has become a major goal. In contrast with flow based works such as those in [START_REF] Christensen | A deformable neuroanatomy textbook based on viscous fluid mechanics[END_REF][START_REF] Miller | On the metrics and Euler-Lagrange equations of computational anatomy[END_REF][START_REF] Thompson | Growth patterns in the developing brain detected by using continuum-mechanical tensor maps[END_REF], our motivation is as in [1, 21-23, 26, 45, 46, 49]. That is, we aim to compute a map that preserves certain pre-defined geometric characteristics of the surfaces. While the literature has mainly attempted to preserve angles (e.g., [START_REF] Angenent | Conformal geometry and brain flattening[END_REF][START_REF] Gu | Genus zero surface conformal mapping and its application to brain surface mapping[END_REF]), and areas (e.g., [START_REF] Chung | Deformation-based surface morphometry applied to gray matter deformation[END_REF][START_REF] Fischl | High-resolution inter-subject averaging and a coordinate system for the cortical surface[END_REF]), we work with geodesic distances (see also [START_REF] Schwartz | A numerical solution to the generalized mapmaker's problem: Flattening nonconvex polyhedral surfaces[END_REF]). Our work is inspired by the mathematically rich literature on Lipschitz minimizing maps and in its connection to the infinite Laplacian. The motivation for using these frameworks will be presented after some brief mathematical introduction below.

Key contributions

Let us conclude this section by explicitly formulating key contributions that make our work stand out from the rich literature in the area of surface, and in particular brain, warping:

-We introduce the theoretical and computational framework of Lipschitz minimizing maps into the area of nonrigid surface and brain mapping. This is the only map, following basic axioms [START_REF] Caselles | An axiomatic approach to image interpolation[END_REF][START_REF] Caselles | An axiomatic approach to scalar data interpolation on surfaces[END_REF], that can simultaneously handle diverse boundary conditions, for example, a mixture of points and curves (see also [START_REF] Glaunes | Diffeomorphic Matching of Distributions: A New Approach for Unlabelled Point-Sets and Sub-Manifolds Matching[END_REF]). -We propose to preserve geodesic distances, being this the first time this is explicitly done for surface warping without intermediate projections. 5 Prior work has concentrated on preserving areas, angles, or Euclidean distances. It is indicated by our preliminary results shown in this paper, the vast literature on geometric surface matching, and the general knowledge in the brain imaging community, that geodesics and intrinsic brain measurements are critical. -All the computations are done intrinsically on the surface, avoiding the distortions commonly introduced by arbitrary parametrizations or projections onto the plane or sphere.

Although the framework as here described is limited to matching natural geodesics distances, additional feature-matching terms could be easily added for example in the form of extra energy terms in the overall cost to be minimized, for which minimizations techniques different than the efficient one here suggested might need to be developed. Also, see for example concluding remarks, additional matching conditions can be simply considered as metric changes, thereby being included in the proposed framework and in the here suggested numerical implementation. To keep the exposition simpler, we do not include these terms in the current presentation.

Formal statement of the problem

Let B 1 and B 2 be two cortical surfaces (2D surfaces in the three dimensional Euclidean space) which we consider smooth and endowed with the metric inherited from IR 3 , so that d B 1 and d B 2 are the geodesic distances measured on B 1 and B 2 , respectively. 6 Let Γ 1 ⊂ B 1 and Γ 2 ⊂ B 2 be subsets which represent features for which a correspondence is already known (these sets could be empty). In general, the sets Γ i are the union of smooth curves traced on the surfaces, e.g., sulcal beds lying between gyri, and/or a union of isolated points. A set of anatomical landmarks that occur consistently in all subjects can be reliably identified using standardized anatomical protocols or automated sulcal labeling techniques (see for example Brain VISA [4] by Mangin and Riviere and SEAL [40] by Le Goualher). The overall goal is to extend the given map between Γ 1 ⊂ B 1 and Γ 2 ⊂ B 2 to the whole surface B 1 and B 2 , minimizing as much as possible the (geodesic) distortion.

Functional anatomy also varies with respect to sulcal landmarks, but sulci typically lie at the interfaces of functionally different cortical regions so aligning them improves the registration of functionally homologous areas. As commonly done in brain warping [START_REF] Toga | Brain Warping[END_REF], we assume that a correspondence between Γ 1 and Γ 2 is prespecified to the map (boundary conditions of the map). In this correspondence, internal point correspondences may be allowed to relax along landmark curves in the final mappings, e.g., [START_REF] Leow | Brain structural mapping using a novel hybrid implicit/explicit framework based on the level-set method[END_REF]. 7To fix ideas let's assume that Γ 1 = ∪ N k=1 x i and Γ 2 = ∪ N k=1 y i , and that the correspondence is given by x i → y i for 1 ≤ i ≤ N . The x i and y i are elements on the surface, such as points or curves. We want to find a (at least continuous) map φ : B 1 → B 2 such that φ(x i ) = y i for 1 ≤ i ≤ N and such that φ produces minimal distortion according to some functional J. One possible way of interpreting this problem is that we are trying to extrapolate or extend the correspondence from Γ 1 to the whole of B 1 in such a way that we achieve small distortion.

A common way to measure the distortion produced by a map φ is by computing the functionals (1 ≤ p < ∞)

J p (φ) = 1 µ(B 1 ) B1 D B1 φ p 2 µ(dx) 1/p (1)
where D B 1 denotes differentiation intrinsic to the surface B 1 and µ is the area measure on B 1 . One immediate idea is then to consider, for a fixed p ∈ (1, ∞), the following variational problem:

Problem 1 (minimize J p ). Find φ ∈ S such that J p (φ) = inf ψ∈S J p (ψ)
, where S is a certain smoothness class of maps φ from B 1 to B 2 such that they respect the given boundary conditions φ(x i ) = y i for all

x i ∈ Γ 1 .
The case p = 2 corresponds to the Dirichlet functional and has connections with the theory of (standard) Harmonic Maps. In more generality, it is customary to call the solutions to Problem 1 p-Harmonic Maps, see for example [START_REF] Eells | A report on harmonic maps[END_REF][START_REF] Fardoun | Heat flow for p-harmonic maps with small initial data[END_REF][START_REF] Jost | Riemannian Geometry and Geometric Analysis[END_REF]. It is easy to show, under mild regularity assumptions, that for a fixed φ, J p (φ) is nondecreasing as a function of p, and that [START_REF] Gilbarg | Elliptic Partial Differential Equations of Second Order[END_REF]

J ∞ (φ) := lim p↑∞ J p (φ) = essup x∈B 1 D B 1 φ(x) 2 , (2)
which is the Lipschitz constant of φ.

In this paper we propose to use the functional J ∞ as a measure of distortion for maps between surfaces and to solve the associated variational problem in order to find a candidate mapping between the surfaces (constrained by the provided boundary conditions). Once again, although our working example are brain surfaces, the here proposed framework is general for non-rigid surface matching.

Let L denote the space of all Lipschitz continuous maps ψ : B 1 → B 2 such that ψ(x i ) = y i for 1 ≤ i ≤ N . We then propose to solve the following problem:

Problem 2 (minimize J ∞ ). Find φ ∈ L such that J ∞ (φ) = inf ψ∈L J ∞ (ψ).
We now argue in favor of this functional.

Why use J ∞ ?

Our first argument is that J ∞ measures distortion in a more global way than any of the J p for p ∈ (1, ∞), since instead of computing an averaged integral quantity, we are looking at the supremum of the local distortions, D B 1 φ(x) 2 . Note also that J p is bounded above by J ∞ under mild regularity assumptions for all p ∈ (1, ∞).

Another element to consider is that this problem is well posed for the kind of general boundary data we want to respect, provided both at curves and isolated points on the surfaces. At least for the case p ≤ 2, this is not true in general, see [START_REF] Caselles | An axiomatic approach to image interpolation[END_REF].

We are then looking for a Lipschitz extension of the map given at Γ 1 whose Lipschitz constant is as small as possible. Let

L(Γ 1 , Γ 2 ) := max x i ,x j ∈Γ 1 d B 2 (yi,yj ) d B 1 (x i ,x j )
, that is, the Lipschitz constant of the boundary data. In general, we have inf ψ∈L J ∞ ≥ L(Γ 1 , Γ 2 ) (we cannot get lower distortion than the one already introduced by the pre-specified boundary conditions). This is related to Kirszbraun's Theorem, which in one of its many guises states that a Lipschitz map f : S → IR D , S ⊂ IR d , has an extension f : IR d → IR D with the same Lipschitz constant as f , see [START_REF] Federer | Geometric Measure Theory[END_REF]. In the same vein, one has Whitney and McShane extensions which apply to the case when the domain is any metric space X and the target is IR. These extensions provide functions that agree with f where boundary conditions are given and preserve the Lipschitz constant throughout X, see for example [START_REF] Aronsson | A tour of the theory of absolutely minimizing functions[END_REF][START_REF] Juutinen | Absolutely minimizing Lipschitz extensions on a metric space[END_REF]. The more general problem of extending f : S → Y (S ⊂ X, X and Y any metric spaces) to all X with the same Lipschitz constant is not so well understood and only partial results are known, see for example [START_REF] Lang | Extendibility of large-scale Lipschitz maps[END_REF][START_REF] Lang | Extensions of Lipschitz maps into Hadamard spaces[END_REF][START_REF] Lang | Kirszbraun's theorem and metric spaces of bounded curvature[END_REF].

The idea then is to keep the distortion at the same order as that of the provided boundary conditions. In general there might be many solutions for the Problem (2). One particular class of minimizers which has recently received a lot of attention is that of absolute minimizers, or absolutely minimizing Lipschitz extensions (AMLE). Roughly speaking, the idea here is to single out those solutions of Problem (2) that also possess minimal local Lipschitz constant, again, see [START_REF] Aronsson | A tour of the theory of absolutely minimizing functions[END_REF] for a general exposition, and [START_REF] Juutinen | Absolutely minimizing Lipschitz extensions on a metric space[END_REF] for a treatment of the case when the domain is any reasonable metric space and the target is the real line. This is the direction we pursue with the computational approach detailed in the next section. 83 Proposed computational approach If we take for example the case of p-Harmonic maps, one way of dealing with the computation of the optimal map φ p is by implementing the geometric p-heat flow associated with the Euler-Lagrange equation of the functional J p , starting from a certain initial condition. As was explained in [START_REF] Mémoli | Solving variational problems and partial differential equations mapping into general target manifolds[END_REF], using an implicit representation for both B 1 and B 2 , we could obtain the partial differential equation PDE p we need to solve in order to find φ p . By taking the formal limit as p ↑ ∞ we would find PDE ∞ , the PDE that characterizes the solution φ ∞ of the (variational) Problem (2). 9 All of this might work if we had a notion of solution for the resulting PDEs. Whereas this is feasible in the case of PDE p for 1 < p < ∞, to the best of our knowledge, there is no such notion of a solution for PDE ∞ . One could of course still persist and try to solve these equations without the necessary theoretical foundations and call these plausible solutions ∞-Harmonic Maps. Nonetheless, this is certainly an interesting line of research.

A different direction is considered in this work. As a guiding example, we first concentrate on the case where B 1 is any closed smooth manifold and B 2 is replaced by IR, as considered in [START_REF] Caselles | An axiomatic approach to scalar data interpolation on surfaces[END_REF] (for scalar data interpolation on surfaces), and in [START_REF] Oberman | Convergent difference schemes for the infinity Laplacian: Construction of absolutely minimizing Lipschitz extensions[END_REF]. In [START_REF] Caselles | An axiomatic approach to scalar data interpolation on surfaces[END_REF], the authors propose to follow a similar path to the one we have just described, and they do not obtain a convergent numerical discretization for the resulting PDE. Meanwhile, in [START_REF] Oberman | Convergent difference schemes for the infinity Laplacian: Construction of absolutely minimizing Lipschitz extensions[END_REF], the author proposes a convergent discretization of the PDE, basing his construction on the original variational problem. We choose to follow this idea as our guiding principle, and extend it to the case where both B 1 and B 2 are surfaces.

We now explain this alternative approach. The basic idea is simple, instead of first obtaining the Euler-Lagrange equations for the energy J ∞ and then discretizing them, we first discretize the energy J ∞ and then proceed to solve the resulting discrete problem. Consider that the domain B 1 is given discretely as a set of (different) points B 1 = {x 1 , . . . , x m } together with a neighborhood relation (i.e., a graph). To fix ideas let's assume the neighborhood relation is a k-nearest neighbors one. Denote, for each 1 ≤ i ≤ m, by N i = {x j 1 , . . . , x j k } ∈ B 1 the set of k neighbors of the point x i . We consider the discrete local Lipschitz constant of the map φ at x i :

L i (φ) := max x j ∈N i d B2 (φ(x i ), φ(x j )) d B1 (x i , x j ) (3)
Upon noting that L i (φ), which measures the local geodesic deformation for the N i neighborhood of the point x i , serves as a discrete approximation to D B 1 φ(x i ) 2 , we see that a possible discretization of the functional J ∞ (φ) is given by the discrete global Lipschitz constant of φ given by max 1≤i≤m L i (φ). The author of [START_REF] Oberman | Convergent difference schemes for the infinity Laplacian: Construction of absolutely minimizing Lipschitz extensions[END_REF] proposed, in the case when B 2 is replaced by IR, solving the discrete version of Problem (2) by the following iterative procedure (here extended for B 2 a surface as in our problem):

-Let φ 0 be an initial guess of the map (see below for details on this for our case).

-For each n ≥ 1, if x i / ∈ Γ 1 , let φ n (x i ) = arg min y∈B2 max j∈N i d B 2 (y, φ n-1 (x j )) d B 1 (x i , x j ) (4) -φ n (x i ) = y i for all n ≥ 0 for x i ∈ Γ 1 .
With computational efficiency related modifications described below, this is the approach we follow in general. The intuition behind this iterative procedure is that, at each point of the domain, we are changing the value of the map in order to minimize the local Lipschitz constant, that is, the local (geodesic) distortion produced by the map. This is in agreement with the notion of AMLEs explained in §2.1. We should remark that since we are using intrinsic distances for the matching, we can let L i (φ) play the role of (the norm of) the displacement field for analyzing the deformation,10 see §5 ahead.

Implementation details

In addition to discretizing the domain B 1 , we also use a discretization B 2 = {y 1 , . . . , y m } of the target space B 2 for our implementation. We endow B 2 with a neighborhood relation given by the k-nearest neighbors of each point. For computational efficiency, we work at all times with two different scales in the discrete domain B 1 . For large data sets, additional scales can be used. We choose a subset F 1 of B 1 such that #F 1 m but still F 1 is an efficient (well separated) covering of B 1 with small covering radius. We do this by using the well known (geodesic) Farthest Point Sampling (FPS) procedure, see [START_REF] Mmoli | A theoretical and computational framework for isometry invariant recognition of point cloud data[END_REF][START_REF] Moenning | Fast marching farthest point sampling[END_REF], which can be efficiently constructed based on optimal computational techniques. Roughly speaking, we apply the iterative procedure on this subset of points only and then extend the map to the rest of the points in the domain B 1 . We now show how to obtain a reasonable initial condition φ 0 and then discuss additional details regarding the implementation of the iterative procedure described in the previous section.

Building the initial condition: We compute, for all (xr,xi) . For this step we use the classical Dijkstra's algorithm for approximating the geodesic distances d B 1 and d B 2 since they might be evaluated at faraway points. This is of course run on the graphs obtained from connecting each point to its k-nearest neighbors. The iterative procedure: After φ 0 is computed for all points in the set F 1 , we run the iterative procedure from the previous section on this set of points. The main modification here is that whereas we still use Dijkstra's algorithm for approximating d B 2 in the target surface, since in the domain we must compute d B1 only for neighboring points (F 1 was chosen to be dense enough), for computational efficiency we can approximate d B1 (x i , x j )

x r ∈ F 1 \Γ 1 , φ 0 (x r ) = arg min y∈B2 max xi∈Γ1 d B 2 (y,y i ) d B 1
x i -x j for x j ∈ N i . We should also point out that for points in F 1 , the neighborhood relation is defined to be that of k-nearest neighbors with respect to the metric on B 1 defined by the adjacency matrix of B 1 . Let φ * : F 1 → B 2 denote the map obtained as the output of this stage. Extension to the whole domain: After we have iterated over points in F 1 until convergence, we extend the map φ * to all points x) . For this step, and since we have already obtained the map for a relatively dense subset, we approximate both d B1 and d B2 by the Euclidean distance. Once again, the motivation for this is just computational efficiency.

x i in B 1 \{F 1 ∪Γ 1 }. This is done by computing φ * (x i ) = arg min y∈B 2 max x∈F 1 ∪Γ 1 d B 2 (y,φ * (x)) d B 1 (x i ,

Examples

In this section we present some computational examples of the ideas presented in previous sections. First, in Figure 1, top, the domain B 1 is a cube (m = 10086) and the target B 2 is a sphere (m = 17982). For the purposes of visualizing the map, we assigned the clown texture (which can be thought of as a function I : B 2 → IR) to the sphere, which can be seen on the bottom-right corner of the figure. The sphere and the cube were concentric and of approximately the same size. We selected F 1 on the cube consisting of 1000 well separated points using the FPS procedure alluded to in §4. Also, we set k = 6 (number of neighbors). We then chose Γ 1 to be the first 100 points of the set and then projected them onto the sphere, obtaining in this way, the corresponding set Γ 2 to use as boundary conditions. We then followed the computational procedure detailed before. The top-left figure shows the composition I • φ * : cube → IR as a texture on the cube. Finally, the top-right and the bottom-left images show the histogram of L i (φ * ) and its spatial distribution in the domain (we paint the cube at each point x i with the color corresponding to L i (φ * )), respectively. Ideally, we would like to obtain a δ-type histogram, meaning that the distances have been constantly scaled. Of course, this is not possible (unless one of the surfaces is isometric to a scaled version of the other), and we attempt to obtain histograms as concentrated as possible. This is quite nicely obtained for this and the additional examples in this paper.

Figure 1, bottom left and right, shows the construction of a map from the unit sphere S 2 into a cortical hemisphere B (B). The boundary conditions consisted of 6 pairs of points. We first took the following 6 points on the sphere Γ 1 = {(±1, 0, 0), (0, ±1, 0), (0, 0, ±1)}. We then constructed the intrinsic distance matrix [d S 2 (p i , p j )] for p i , p j ∈ Γ 1 . Finally, we chose 6 points {q 1 , . . . , q 6 } = Γ 2 in B such that max i =j

d B (qi,qj ) d S 2 (p i ,p j )
was as close as possible to 1 π diam(B). We painted B with a texture I H depending on its mean curvature so as to more easily visualize the sulci/crests: If H(x) stands for mean curvature of B at x, then I H (x) = (H(x) -min x H(x))

2 . See the caption for more details.

The example in Figure 2 is about computing a map Φ from a subject's left hemisphere B 1 to another subject's left hemisphere B 2 . The boundary conditions were constructed in a way similar to the one used for the previous example, but in this case, 300 points were chosen. Note that, if available, hand traced curves could be used as commonly done in the literature (in other words, more anatomical/functional oriented boundary conditions). In the first two rows we show 4 different views of each cortical surface, and in the third row we show B 1 colored with the values of L i (Φ) which we interpret as a measure of the local deformation of the map needed to match B 1 to B 2 . See the caption for more details.

Concluding remarks

In this paper we have introduced the notions of minimizing Lipschitz extensions into the area of surface and brain warping. These maps provide a more global constraint than ordinary p-harmonic ones, and allow for more general boundary conditions. The proposed computational framework leads to an efficient surface-tosurface warping algorithm that avoids distorting intermediate steps that are common in the brain warping literature. We are currently investigating the use of this new warping technique for creating population averages and applying it to disease and growth studies. In earlier work, the Jacobian of a deformation mapping over time has been used to map the profile of brain tissue growth and loss in a subject scanned serially (tensor-based morphometry [START_REF] Chung | Tensor-based brain surface modeling and analysis[END_REF][START_REF] Thompson | Growth patterns in the developing brain detected by using continuum-mechanical tensor maps[END_REF]). The discrete local Lipschitz constants of our computed mappings also provide a useful index of deformation that can be analyzed statistically across subjects. The framework here introduced can also be applied in 3D for volumetric warping and with weighted geodesic distances instead of natural ones, or different metrics in general, to include additional geometric characteristics in the matching. As recently shown in [START_REF] Leordeanu | A spectral technique for correspondence problems using pairwise constraints[END_REF], the use of pairwise distances is of importance of other matching and computer vision tasks beyond the ones discussed in this paper. Results in these directions will be reported elsewhere.

Finally, a further validation and possible different application of the framework here proposed is presented in Figure 3. Here, the 3D shape on the left is mapped onto its given isometric 3D shape on the right. We used 2000 points per shape, with 20 pointwise boundary conditions. The bottom-left shows the shape on the top-left, colored according to the corresponding local Lipschitz constant. The bottom-right shows the histogram of the local Lipschitz constant, centered around one (average is 1.1917). The small deviation from the ideal case for isometric shapes, an histogram with all the points exactly at one, is expected due to the sampling (the shapes were first bended and then sampled), the algorithmic discretization, perturbations in the boundary conditions matches, and the extrapolation from a limited number of points. The detailed exploitation of this framework for metric-invariant shape warping and recognition, with side information in the form of boundary conditions, will be reported elsewhere. On the middle, the texture map corresponding to IH (x) as described in the text is used. On the right, the texture is max(IH (x), δ) for a user selected value of the threshold δ, which highlights the gyral crests.

Fig. 1 .

 1 Fig. 1. Top: Artificial example of the proposed warping algorithm. From top to bottom and left to right: The domain surface, with a picture painted on it to help in visualizing the computed map; histogram of the Lipschitz constant (note how it is concentrated around a single value); color coded distribution of the Lipschitz constant for the computed map; and mapped texture following the computed map. Bottom, left and right: Example of mapping between the cortex and a sphere. The order is the same as in the previous figure, but now the domain and target surface are colored with a curvature-based color code. Note once again the concentration of the Lipschitz constant for the computed map.On the middle, the texture map corresponding to IH (x) as described in the text is used. On the right, the texture is max(IH (x), δ) for a user selected value of the threshold δ, which highlights the gyral crests.

Fig. 2 .

 2 Fig. 2. Warping between the cortical surfaces of two brains. In the first row, we show 4 views of B 1 : posterior, medial, lateral and directly viewing the occipital cortex. The corresponding 4 views of B 2 are shown on second row. In the third row, we show B 1 with texture I(x i ) = L i (Φ) which can interpreted as a measure of local deformation needed to match x i ∈ B 1 to Φ(x i ) ∈ B 2 . Relatively little deformation (blue colors) is required to match features across subjects on the flat interhemispheric surface (second image in the second row). This is consistent with the lower variability of the gyral pattern in the cingulate and medial frontal cortices. By contrast, there is significant expansion required to match the posterior occipital cortices of these two subjects, especially in the occipital poles which are the target of many functional imaging studies of vision. The final panel in the figure shows the corresponding histogram for L i (φ), the local Lipschitz constants of the map.

Fig. 3 .

 3 Fig.3. Mapping results for isometric invariant 3D shapes.

This includes groups at JHU, UCLA, U. Penn., INRIA, MGH, GATECH, Harvard-BW, and the University of Florida, to name just a few.

After the original version of this paper from March 2005, http://www.ima.umn.edu/preprints/jan2006/jan2006.html, and motivated in part by our work on Gromov-Hausdorff for isometric shape comparison,[START_REF] Mmoli | A theoretical and computational framework for isometry invariant recognition of point cloud data[END_REF], and by their work on multidimensional scaling (MDS), the authors of[5,

6] proposed a gradient descent technique to solve the classical STRESS function from MDS, where the attributes they select are geodesic distances.[START_REF] Bronstein | Efficient computation of isometry-invariant distances between surfaces[END_REF] The work developed below applies to any desired intrinsic distances, where we could for example use non-uniform metrics on the surface as in[START_REF] Bartesaghi | A system for the generation of curves on 3D brain images[END_REF][START_REF] Mémoli | Distance functions and geodesics on submanifolds of R d and point clouds[END_REF]. For simplicity of the presentation, we restrict the discussion to the usual geodesic distance.

Extending the framework here introduced, in particular the matching of metric measurements, to handle flexibility in the boundary conditions, e.g.,[START_REF] Glaunes | Diffeomorphic Matching of Distributions: A New Approach for Unlabelled Point-Sets and Sub-Manifolds Matching[END_REF], is an interesting direction to pursue.

As commonly happens in the surface and brain warping literature, the map cannot be guaranteed to be invertible. We have not experienced any problems with this in our experiments.

The case when the domain is a subset of IR d and the target is the real line leads to the so called infinity Laplacian, see[START_REF] Aronsson | A tour of the theory of absolutely minimizing functions[END_REF][START_REF] Crandall | Optimal Lipschitz extensions and the infinity Laplacian[END_REF][START_REF] Jensen | Uniqueness of Lipschitz extensions: Minimizing the sup norm of the gradient[END_REF].

One can imagine a situation in which two isometric surfaces are matched by our algorithm such that Li(φ) = 1 for all i, but the displacement field x i -φ(x i ) is large since there may be no rigid motion that aligns the two surfaces. One simple example is a flat sheet of paper and the same sheet slightly bent.

Mathematical Foundations of Computational Anatomy (MFCA'06)

This work is partially supported by DARPA, the Office of Naval Research, the National Science Foundation, the National Institutes of Health, and the National Geospatial-Intelligence Agency.