Maxime Leroux
email: mleroux@anl.gov

Ion Errea

Mathieu Le Tacon

Sofia-Michaela Souliou

Gaston Garbarino

Laurent Cario

Alexey Bosak

Francesco Mauri

Matteo Calandra
email: matteo.calandra@upmc.fr

Pierre Rodière
email: pierre.rodiere@neel.cnrs.fr

Strong anharmonicity induces quantum melting of charge density wave in 2H-NbSe2 under pressure

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

The pressure and temperature dependence of the phonon dispersion of 2H-NbSe2 is measured by inelastic X-ray scattering. A strong temperature dependent soft phonon mode, reminiscent of the Charge Density Wave (CDW), is found to persist up to a pressure as high as 16 GPa, far above the critical pressure at which the CDW disappears at 0 K. By using ab initio calculations beyond the harmonic approximation, we obtain an accurate, quantitative, description of the (P,T) dependence of the phonon spectrum. Our results show that the rapid destruction of the CDW under pressure is related to the zero mode vibrations -or quantum fluctuations -of the lattice renormalized by the anharmonic part of the lattice potential. The calculations also show that the low-energy longitudinal acoustic mode that drives the CDW transition barely contributes to superconductivity, explaining the insensitivity of the superconducting critical temperature to the CDW transition.

The interplay between charge density wave order, i.e. a static modulation of the electronic density close to the Fermi level, and superconductivity has attracted much attention 1,[START_REF] Le Tacon | Inelastic X-ray scattering in YBa2Cu3O6.6 reveals giant phonon anomalies and elastic central peak due to charge-density-wave formation[END_REF] . This is the central issue of a long standing debate in simple transition metal dichalcogenides without strong electronic correlations, such as 2H-NbSe 2 [START_REF] Moncton | Neutron scattering study of the charge-density wave transitions in 2H-NbSe2 and 2H-TaSe2 by Neutron Scattering[END_REF][START_REF] Wilson | Charge-density waves and superlattices in the metallic layered transition metal dichalcogenides[END_REF] . At T CDW =33.5 K, 2H-NbSe 2 undergoes a second order phase transition towards an incommensurate CDW phase which coexists with superconductivity below T c = 7.2 K. The anisotropy and temperature dependence of the electronic band structure has been widely studied by ARPES and STM measurements [5][6][7][START_REF] Guillamon | Intrinsic atomic-scale modulations of the superconducting gap of 2H-NbSe2[END_REF][9][START_REF] Rahn | Gaps and kinks in the electronic structure of the superconductor 2H-NbSe2 from angle-resolved photoemission at 1 K[END_REF][START_REF] Arguello | Quasiparticle Interference Quasiparticle Interactions and the Origin of the Charge Density Wave in 2H-NbSe2[END_REF] . The CDW is coupled to a periodic lattice distortion through a strong electronphonon coupling. The transition is associated with a softening of a longitudinal acoustic phonon mode as temperature is lowered to T CDW [START_REF] Weber | Extended Phonon Collapse and the Origin of the Charge-Density Wave in 2H-NbSe2[END_REF] . It was already noticed that the high order phonon fluctuations and strong electronphonon interactions explain some of the key features of the formation of the CDW in this system [START_REF] Inglesfield | Bonding and phase transitions in transition metal dichalcogenide layer compounds[END_REF][START_REF] Varma | Strong-Coupling Theory of Charge-Density-Wave Transitions[END_REF] . In the specific case of NbSe 2 at ambient pressure, the reduction of the phonon lifetime through the electron-phonon coupling is a key ingredient on the temperature dependence of the phonon dispersion [START_REF] Flicker | Charge order from orbitaldependent coupling evidenced by NbSe2[END_REF] .

At ambient pressure, published ab initio calculations reproduce successfully the lattice instability on a wide part of the Brillouin zone around the experimentally observed q CDW [START_REF] Weber | Extended Phonon Collapse and the Origin of the Charge-Density Wave in 2H-NbSe2[END_REF][START_REF] Calandra | Effect of dimensionality on the charge-density wave in few-layer 2H-NbSe2[END_REF][START_REF] Weber | Optical phonons and the soft mode in 2H-NbSe2[END_REF] . However, these calculations are carried out at the harmonic level, and therefore at T=0K. The temperature dependence can be qualitatively reproduced using a Gaussian smearing of the Fermi surface which reduces the contribution of the electron-phonon interaction. Quantitatively however, the temperatures that would effectively yield this smearing are unphysical, as they exceed the experimentally observed one by several orders of magnitude.

The quantitative failure of the harmonic ab-initio calculations is even worse under pressure. Experimentally, by applying an hydrostatic pressure, T CDW is driven to 0 K and the CDW instability disappears above the critical pressure P CDW (T = 3.5 K)=4.6 GPa 18 . Around this quantum phase transition superconductivity remains essentially unaffected, as shown on the (T,P) phase diagram in Fig. 1 [START_REF] Berthier | Evidence for a connection between charge density waves and pressure enhancement of superconductivity in 2H-NbSe2 Solid State Communications[END_REF][START_REF] Suderow | Pressure induced effects on the fermi surface of superconducting 2H-NbSe2[END_REF] . However, harmonic ab initio calculations fail to describe the ground state of NbSe 2 as a function of pressure (see SI and 21 ). Indeed, they predict (see Fig. 3 and SI) a CDW instability up to P C = 14 GPa, a pressure that largely exceeds the experimentally observed one. This persistence of the CDW instability in the harmonic calculations is reminiscent of the effect of iso-valent and iso-electronic substitution of Se with lighter S. The harmonic calculation also predicts a CDW instability for NbS 2 that is not observed experimentally [START_REF] Leroux | Anharmonic suppression of charge density waves in 2H-NbS2[END_REF] .

Taken together, these facts indicate that for an accurate description of the ground state (even at T=0K) of these systems, calculations beyond the standard harmonic approximation must be carried out. The impact of anharmonicity has recently been highlighted in different systems. To mention only a few, the thermodynamical properties of the metal-insulator transition in strongly correlated VO 2 [START_REF] Budai | Metallization of Vanadium dioxide driven by large phonon entropy[END_REF] , the transport properties of thermoelectric PbTe [START_REF] Delaire | Giant anharmonic phonon scattering in PbTe[END_REF] or the electron-phonon interaction above the Verwey transition in magnetite Fe 3 O 4 [START_REF] Hoesch | Anharmonicity due to Electron-Phonon Coupling in Magnetite[END_REF] are all strongly influenced by anharmonicity. Even at low temperatures, where its effects are expected to be weaker, it must be taken into account to explain the pressure dependence of the crystal structure of Calcium 26 , the inverse isotope effect in superconducting PdH 27 , or the phonon spectra at the proximity of the liquid-solid transition of Helium at 0K [START_REF] Glyde | Excitations in liquid and solid Helium[END_REF] . Finally, anharmonicity is crucial for understanding the ground state of a system when a quantum electronic phase transition coupled to the lattice occurs, as in e.g. ferroelectrics [START_REF] Rowley | Ferroelectric quantum criticalityNat[END_REF] or CDW compounds [START_REF] Weber | Extended Phonon Collapse and the Origin of the Charge-Density Wave in 2H-NbSe2[END_REF][START_REF] Varma | Strong-Coupling Theory of Charge-Density-Wave Transitions[END_REF][START_REF] Leroux | Anharmonic suppression of charge density waves in 2H-NbS2[END_REF][START_REF] Grüner | Density Waves in Solids[END_REF] . In most cases however, these effects are computationally too expensive to be evaluated in standard first-principle calculations, and therefore often disregarded.

Here, we report the investigation of the temperature and pressure dependence of the phonon spectra of NbSe 2 by using high resolution inelastic X-ray scattering (IXS) from a crystal in a diamond anvil cell, up to pressures as large as 16 GPa at T = 20 K. We observe that the strong temperature dependence of the soft phonon mode is reduced by pressure, but still present even at 16 GPa.

We show that the temperature and pressure dependence of the phonon spectra can be accurately accounted for when the effects of anharmonicity are explicitly included in the calculation. Anharmonic effects were calculated within the newly developed stochastic selfconsistent harmonic approximation (SSCHA), yielding an unprecedented agreement with the experimentally determined temperature dependence of the phonon spectra at various pressures. We show that at low temperatures, for pressures between P CDW and P C , the CDW is de-stroyed by quantum fluctuations, which easily overcome the double-well lattice potential and, thus, are strongly affected by the anharmonic part of the potential. Surprisingly, we observe that anharmonic effects dominate the low-energy phonon spectrum of 2H-NbSe 2 in a large pressure range, extending at least up to 16 GPa >> P CDW and down to the lowest temperatures. Finally, we demonstrate that a large electron-phonon interaction, mostly due to optical modes rather than to the soft longitudinal acoustic mode which drives the CDW instability, contributes substantially to superconductivity in 2H-NbSe 2 , and naturally explains its insensitivity to the occurrence of a CDW.

Single crystals of 2H-NbSe 2 were grown using the vapor growth technique in sealed quartz tubes with iodine as a transport agent 3 . Typical dimensions of the crystals are 100×100×50 µm 3 ( a× b× c). Pressure was generated using a diamond anvil cell, loaded with helium as a pressure transmitting medium to ensure highly hydrostatic conditions. Temperature was lowered using a customdesigned [START_REF] Wilson | Charge-density waves and superlattices in the metallic layered transition metal dichalcogenides[END_REF] He cryostat, and pressure was varied in-situ at low temperature using a helium-pressurized membrane. Each cell contained two rubies to monitor the pressure insitu using the fluorescence technique. Measurements of the phonon dispersion of 2H-NbSe 2 were taken between room temperature and 20 K, and for pressures ranging from 0.2 to 16 GPa.

IXS measurements were carried out on beamline ID28 at the ESRF. The x-ray beam was aligned along the caxis of the crystal. We used the (9,9,9) reflection on the backscattering monochromator with an incident energy of 17.794 keV, and a corresponding energy resolution of 1.3 meV HWHM (lorentzian fit of the elastic peak). The incident beam was focused by a multilayer mirror into a spot of 100×60µm (width × height). We used 20×60mm (width × height) slits at a 7 m distance giving a momentum resolution of 0.014 Å-1 (a * × c * ≈ 2.1 × 0.50 Å-1 ) in the (H,0,L) plane, and 0.042 Å-1 (b * ×sin(60 • ) ≈ 1.8 Å-1 ) in the perpendicular direction. We measured the phonon dispersion along the ΓM direction close to Γ 200 (i. e.

(2-h, 0, 0)), where both longitudinal optical and acoustic like soft phonon branches can be observed. According to the balance of the dynamical structure factor of the 2 phonon modes and so of the inelastic scattering function S(q,ω), the observation of one of the two branches can be selected. For each transfer wave vector investigated, an energy scan up to 15meV was performed.

All the spectra were fitted using the standard Levenberg-Marquadt algorithm, with the following free parameters : position and amplitude of the elastic peak; position, linewidth and amplitude of the phonon; and a constant background. We make the standard assumptions [START_REF] Fåk | Phonon line shapes and excitation energies[END_REF][START_REF] Burkel | Phonon spectroscopy by inelastic x-ray scattering[END_REF] that the elastic peak is Lorentzian, and that the phonon lineshape can be modeled by a damped harmonic oscillator convoluted by a normalized Lorentzian with a width corresponding to the resolution of the detector.

Ab initio calculations were performed within density functional perturbation theory and the general-ized gradient approximation 33 making use of the Quantum ESPRESSO [START_REF] Giannozzi | QUANTUM ESPRESSO: a modular and open-source software project for quantum simulations of materials[END_REF] code. An ultrasoft pseudopotential (norm-conserving) for Nb (Se), a 35 Ry energy cutoff, and a 24×24×8 mesh for the electronic integrations were used. An Hermitian-Gaussian smearing of 0.01 Ry was used. Harmonic dynamical matrices were calculated within linear response in a grid of 6×6×4 q points. The small smearing parameter and the dense electronic mesh used are needed to converge the phonon frequency of the longitudinal acoustic mode at the CDW wave vector (cf. SI). The stochastic self-consistent harmonic approximation (SSCHA) [START_REF] Errea | First-Principles Theory of Anharmonicity and the Inverse Isotope Effect in Superconducting Palladium-Hydride Compounds[END_REF][START_REF] Errea | Anharmonic free energies and phonon dispersions from the stochastic self-consistent harmonic approximation: Application to platinum and palladium hydrides[END_REF] calculations were performed using a 3×3×1 supercell, yielding anharmonic dynamical matrices in the commensurate 3×3×1 grid. The difference between the SSCHA dynamical matrices and the harmonic dynamical matrices in the 3×3×1 grid was Fourier interpolated to the points of the finer 6×6×4 grid. Adding the harmonic dynamical matrices in the 6×6×4 grid to the result of the interpolation, the SSCHA dynamical matrices in the finer 6×6×4 mesh were obtained. The experimental lattice parameters [START_REF] Feng | Order parameter fluctuations at a buried quantum critical point[END_REF] with relaxed internal coordinates were used because in these conditions the results weakly depend on the exchange-correlation potential chosen.

The electron-phonon calculations were performed by using maximally localized Wannier functions [START_REF] Marzari | Maximally localized generalized Wannier functions for composite energy bands[END_REF][START_REF] Souza | Maximally localized Wannier functions for entangled energy bands[END_REF] and interpolation of the electron-phonon matrix elements as in Ref. [START_REF] Calandra | Adiabatic and nonadiabatic phonon dispersion in a Wannier function approach[END_REF] . We used 14 Wannier functions (d-states of Nb and Se p z states) and a 6 × 6 × 4 electron-momentum grid. We used the Fourier interpolated SSCHA dynamical matrices in the electron-phonon coupling calculations. The average electron-phonon coupling was interpolated using electron and phonon momentum grids of 24 × 24 × 8 randomly displaced from the origin.

The temperature dependence of the low energy phonon dispersions recorded at pressures of 3.5 and 16 GPa are shown in Fig. 2 (some of the corresponding raw data is shown in the supplementary materials). At T = 20 K and P = 3.5 GPa, the system is very close to the quantum critical regime reported in ref. 18. No particular enhancement of the elastic line intensity close to the CDW wavevector q CDW ≈ (1.67, 0, 0) is observed. On the other hand, a colossal temperature dependent softening of the longitudinal acoustic phonon mode is seen : from ∼ 9 to ∼ 1.5 meV between room temperature and 20 K. While this mode is expected to condense around ∼ 12 K 18 , this temperature could not be reached at this pressure with the current experimental set-up. As we increase pressure to 16 GPa, the amplitude of the phonon softening is significantly reduced, but remains sizeable as the mode's energy almost decreases by a factor of 2 from 300 to 20 K. So far, most of the theoretical and experimental work on NbSe 2 has focused on the interplay between the electronic susceptibility and the electron-phonon interaction as possible mechanisms for the CDW formation [START_REF] Johannes | Fermisurface nesting and the origin of the charge-density wave in NbSe2[END_REF][START_REF] Zhu | Classification of charge density waves based on their nature Proc[END_REF] . The large temperature dependence indicates that strong anharmonic effects survive over a pressure range that is way larger than the experimentally reported CDW sta- bility range (0 < P < P CDW ) and might therefore play an important role that was up to now largely disregarded. In Fig. 3, we plot the square of the soft phonon frequency measured at T = 20 K at q CDW as a function of pressure. The linear dependence indicates that the mode hardens as √ P -P CDW . We also plot the square of the phonon frequency calculated within the harmonic approximation, which is negative -i.e. indicating that the system is unstable against CDW formation -up to 14 GPa >> P CDW . This is reminiscent of the situation encountered at ambient pressure in NbS 2 , where we have observed a large phonon softening upon cooling, insufficient however to induce the CDW formation despite the predictions of the harmonic ab-initio calculation [START_REF] Leroux | Anharmonic suppression of charge density waves in 2H-NbS2[END_REF] .

To quantify the effect of anharmonicity and check whether it is responsible for the destabilization of the CDW between P CDW and 14 GPa, we have carried out a series of calculations for the pressures and temperatures experimentally investigated, within the SSCHA approach [START_REF] Errea | First-Principles Theory of Anharmonicity and the Inverse Isotope Effect in Superconducting Palladium-Hydride Compounds[END_REF][START_REF] Errea | Anharmonic free energies and phonon dispersions from the stochastic self-consistent harmonic approximation: Application to platinum and palladium hydrides[END_REF] (see SI) which allows us to access directly the anharmonic free energy of the system, with full inclusion of the anharmonic potential terms.

The results are superimposed to our experimental dispersions in Fig. 2. Not only is the instability suppressed as experimentally observed, but the agreement between the measured and calculated dispersions over the entire pressure range is remarkable compared to harmonic ab-initio calculation. A small overestimation of phonon calculation is nevertheless observed in the range around q CDW .

A key result here is that, even at low temperature, anharmonic terms play a key role in the destruction of the harmonically stable long-range CDW order. This is surprising at first glance since anharmonic effects usually manifest themselves at finite temperature, through the effect of many-body interactions between thermally excited phonons. These account for the lattice thermal expansion and for the decay of the phonons lifetime at high temperatures, but are usually strongly suppressed at low temperatures. On the other hand, given that atoms In all figures the grey background denotes the integration area used to infer that the soft acoustic longitudinal mode contributes 17% of the total electronphonon coupling.

vibrate even at zero temperature due to the zero-point fluctuations, their motion can in principle be affected by the anharmonic potential. This appears to be the case in this system where, below 14 GPa, the harmonic part of the potential that drives the lattice instability is overwhelmed by the anharmonic part. The balance between these two parts of the potential tells whether there is a lattice instability (driven by the CDW) or not. This bears stark analogy with the destruction of long range magnetic order from spin-spin quantum fluctuations in low dimensional quantum magnets [START_REF] Balents | Spin liquids in frustrated magnets[END_REF] .

Finally we note that these many-body effects are relevant not only at q = q CDW , but extend up to the M point, the boundary of the Brillouin zone at q = a * /2, where two low-energy modes, the longitudinal acoustic mode that drives the CDW transition and an optical mode, are substantially hardened with respect to the result of the harmonic calculation.

The importance of the electron-phonon interaction in the CDW formation has been emphasized [START_REF] Weber | Extended Phonon Collapse and the Origin of the Charge-Density Wave in 2H-NbSe2[END_REF][START_REF] Johannes | Fermisurface nesting and the origin of the charge-density wave in NbSe2[END_REF] , but its role in superconductivity remains debated, as is, more generally, the interplay between these two orders. Indeed, the absence of reliable ab initio phonon calculations in the high symmetry phase, has prevented a calculation of the superconducting critical temperature. Having achieved a detailed understanding of the low energy part of the experimental phonon spectra and having shown how anharmonicity can suppress the long range ordered CDW phase, we can now address the superconducting properties of NbSe 2 by calculating the Eliashberg function and the integrated electron-phonon coupling (see SI) as a function of pressure. As of now, this calculation can only be done in the high pressure phase, where the number of atoms per unit cell is low enough. We show in Fig. 4 the results of the calculation at 7 GPa (data for other pressures are shown in the SI). We find that the average electron-phonon coupling is as large as λ = 1.16 (see Tab. I for other pressures), demonstrating that NbSe 2 is a strong coupling superconductor. Previous works [START_REF] Arguello | Quasiparticle Interference Quasiparticle Interactions and the Origin of the Charge Density Wave in 2H-NbSe2[END_REF][START_REF] Weber | Extended Phonon Collapse and the Origin of the Charge-Density Wave in 2H-NbSe2[END_REF][START_REF] Calandra | Effect of dimensionality on the charge-density wave in few-layer 2H-NbSe2[END_REF][START_REF] Weber | Optical phonons and the soft mode in 2H-NbSe2[END_REF][START_REF] Gor'kov | Strong electron-lattice coupling as the mechanism behind charge density wave transformations in transition-metal dichalcogenides[END_REF] identified the large electron-phonon coupling of the soft acoustic mode at q CDW as the mechanism responsible for the CDW formation. For this reason, it is natural to expect that this soft mode also plays a key role in superconductivity. However, its contribution to the average electron-phonon coupling is less than 17% (as it can be inferred from the integrated λ(ω) in Fig. 4), meaning that it has only a marginal role in superconductivity. Indeed, the strong coupling of this mode is very localized around q CDW , and its contribution to λ averages out when integrating over the Brillouin zone. The modes sustaining superconductivity are given by the two main contributions to the Eliashberg function: a broad peak centred around ≈ 11.5 meV and a second one in the 22 -26 meV region. These two features contribute to ≈ 70% of the total electron-phonon coupling. Further insights are obtained by decomposing the phonon density of states (DOS) in atomic vibrations along different directions. As the DOS and α 2 F (ω) are very similar, the decomposition will also apply to the Eliashberg function. The low energy feature at 11.5 meV is mainly due to Nb and Se vibrations parallel to the Nb plane, corresponding to the flat phonon band along ΓM suffering from a large anharmonic correction, while the features in the 22 -26 meV range are mostly attributed to out-of-plane and in-plane Se displacements with a non negligible Nb component. In both modes sustaining superconductivity there is a substantial Se component, upturning the conventional wisdom that in superconducting dichalcogenides the chalcogene is accessory while the transition metal plays a key role 4 . Our work demonstrates that the superconducting properties strongly depend on both the transition metal and the chalcogene, offering a natural explanation for the different relation between CDW and superconductivity encountered in different dichalcogenides [START_REF] Tissen | Pressure dependence of superconducting critical temperature and upper critical field of 2H-NbS2[END_REF][START_REF] Kusmartseva | Pressure Induced Superconductivity in Pristine 1T-TiSe2[END_REF][START_REF] Sipos | From Mott state to superconductivity in 1T-TaS2[END_REF] . Moreover, as optical phonon modes contribute strongly to superconductivity in 2H-NbSe 2 , while the CDW is determined by the softening of the longitudinal acoustic mode, the superconducting T c is insensitive to the presence of the CDW, explaining the phase diagram in Fig. 1.

Our work points toward the key role of anharmonicity in destroying the CDW order not only through the usual thermal phonon-fluctuations, but also through quantum phonon-fluctuations, and on a very wide pressure range above the quantum phase transition. Our results, in this respect, are fundamental and relevant for a very large set of materials exhibiting a CDW instability -or more generally a second order electronic quantum phase transition coupled to the lattice -such as transition metal dichalcogenides, cuprates, Bechgaard salts and transition metal bronzes.

FIG. 1 .

 1 FIG. 1. Phase diagram of NbSe2. Squares indicate the limit of the CDW phase measured in Ref. 19 with the addition of the point at 4.6 GPa from Ref. 18 , triangles indicate the superconducting Tc measured in 20 , and circles indicate the pressure and temperature of our IXS measurements. Inset : crystallographic structure of 2H-NbSe2 3 .

FIG. 2 .

 2 FIG.2. Phonon dispersion at 3.5 GPa (a) and 16 GPa (b). Closed (resp. open) dots represent the experimentally determined longitudinal acoustic (longitudinal optic) soft phonon dispersions, and lines represent the results of phonon calculations. Close to the Γ point, the phonon dispersion shows clearly an optical branch. Closer to the M point, the energy of the phonon observed correspond to the acoustic like branch calculated. Colored thick and thin lines represent the SSCHA results for the longitudinal acoustic and optical anharmonic phonons, respectively, at several temperatures. The thin black lines represent other phonon modes present in this energy range calculated within the SSCHA at 20 K, which are very harmonic and barely depend on temperature (see SI).

FIG. 3 .

 3 FIG.3. Pressure dependence of the square of the measured energy of the longitudinal acoustic soft phonon at the CDW ordering wavevector and at 20 K (red dots) compared to the square of the energy calculated in the harmonic approximation (black diamonds). Lines are linear fits. The red region indicates the range of pressure where a CDW is observed experimentally at 20 K. The orange region indicates the range of pressure where harmonic calculations predict a CDW ground state.

FIG. 4 .

 4 FIG. 4. (a) Calculated phonon dispersion at 7 GPa and 0 K with the SSCHA. The size of the bars is proportional to the electron-phonon contribution to the phonon linewidth (the bar is twice the FWHM). (b) The Eliashberg function α 2 F (ω) and the integrated electron-phonon coupling λ(ω).(c) The phonon density of states decomposed into the different atoms.In all figures the grey background denotes the integration area used to infer that the soft acoustic longitudinal mode contributes 17% of the total electronphonon coupling.

TABLE I .

 I Calculated λ, logarithmic frequency average ω log , and Tc values for 2H-NbSe2 at 3.5, 7 and 16 GPa. The McMillan equation with µ * = 0.16 is used to calculate Tc. Experimental data are taken from Ref.[START_REF] Suderow | Pressure induced effects on the fermi surface of superconducting 2H-NbSe2[END_REF] .

	P (GPa) λ ω log (meV) Tc calc. (K) Tc exp. (K)
	3.5	1.28	12.6	11.3	7.8
	7	1.16	13.7	10.5	8.2
	16	0.91	16.4	7.8	7.8

ACKNOWLEDGEMENTS

Sabrina Salmon-Bourmand is thanked for her help in the synthesis of 2H-NbSe 2 single crystals. P.R. and M.L. acknowledge financial support from the French National Research Agency through Grant No. ANR-12-JS04-0003-01 SUBRISSYME. M.C, F.M. and I.E. acknowledge support from the Basque Government (Grant No. BFI-2011-65), Spanish Ministry of Economy and Competitiveness (FIS2013-48286-C2-2-P,) Graphene Flagship and ANR (contract ANR-13-IS10-0003-01) and Prace (2014102310). Calculations were performed at IDRIS, CINES, DIPC and at CEA TGCC.


G Goavec-Merou

Y Yakoubi

R Couturier

M Lenczner

J M Friedt

F Yang

FPGA Implementation of Diffusive Realization for a Distributed Control Operator

Keywords: Distributed Control, FPGA, Diffusive Realization, Semidecentralized Control, Fine Grained Parallel Computer Architectures

We focus on the question of real-time computation for optimal distributed filtering or control applicable to MEMS Arrays. We present an algorithm for the realization of a linear operator solution to a functional equation through its application to a Lyapunov operatorial equation associated to the heat equation in one dimension. It is based on the diffusive realization, and turns to be well suited for fined grained parallel computer architecture as Field Programmable Gate Arrays (FPGA). An effective FPGA implementation has been successfully carried out. Here, we report the main implementation steps and the final measured performances.

I. INTRODUCTION

Our concern when developing the method presented in this paper relates to embedded intensive real time computation based on fine grained parallel computer architectures as Field-Programmable Gate Arrays (FPGA). Here we address the general problem of realization of a linear operator u → z = P u in infinite dimensional spaces. Our method is presented and illustrated on the example of the operator P, solution to the Lyapunov equation,

d 2 dx 2 P u + P d 2 dx 2 u = Qu (1)
in ω = (0, 1) for all u vanishing at the boundary of ω, where Q is another linear operator. It is issued from optimal filtering or control theory of the heat equation ∂T ∂t -∂ 2 T ∂x 2 = q in ω. The method is in the same time fast and suitable for implementation on semi-decentralized architectures where communications are mainly between neighbor unit cells. Its formulation relates to an interesting technique, the diffusive representation, applicable to causal operator realizations which is well developed in the context of time operators. One of the main recognized advantages of this approach is its low computational cost, see the papers of G. Montseny et al., e.g. [START_REF] Laudebat | Modeling and optimal identification of pseudodifferential electrical dynamics by means of diffusive representation -Part 1: Modeling[END_REF], and of D. Matignon et al., e.g. [START_REF] Hélie | Criterion design for optimizing low-cost approximations of infinite-dimensional systems: towards efficient real-time simulation[END_REF], for representations of various pseudodifferential operators and for their approximation. Those of C. Lubich et al., e.g. [START_REF] López-Fernández | ast Runge-Kutta approximation of inhomogeneous parabolic equations[END_REF], apply a similar idea to convolution operators and they develop optimized numerical methods.

In this paper, we briefly present the diffusive representation of the solution P to (1) together with algorithms for their approximation implementable in semi-decentralized architectures. We refer the interested reader to the general theory shortly introduced in [START_REF] Lenczner | Diffusive realization of operator solutions of certain operational partial differential equations[END_REF] and detailed in [START_REF] Yakoubi | Deux Méthodes d'Approximation pour le Contrôle Optimal Semi-Décentralisé des Systèmes Distribués[END_REF]. Then, the core of the paper is devoted to FPGA implementation using integer computation. We start by a rescaling of the algorithms guarantying that all intermediary and final results be in a prescribed interval. Then we report the results of our study on the effect of quantization. After these preparation steps, the FPGA mapping is analyzed for three different possible schemes, namely: Iterative, Pipelined and Parallel. The Parallel solution has been implemented in the Spartan3A by Xilinx, and we provide detailed results regarding resource and time consumptions.

II. DIFFUSIVE REPRESENTATION OF P

We consider the kernel formulation of the operator P ,

P u(x) = Z 1 0 p(x, y)u(y) dy,
and the decomposition z = z + +z -into causal and anti-causal parts,

z + = Z x 0 p(x, y)u(y) dy and z -= Z 1 x p(x, y)u(y) dy.
The kernel p is the unique solution to the boundary value problem -∆p = q in the square (0, 1) × (0, 1), and p = 0 on the square boundary (0, 1) × (0, 1),

where q is the kernel of Q. The realization of z + and of z -may be formulated thanks to the diffusive representation, see [START_REF] Lenczner | Diffusive realization of operator solutions of certain operational partial differential equations[END_REF] and [START_REF] Montseny | Représentation diffusive[END_REF], in the form

z + (x) = Z R µ + (x, ξ)ψ + (x, ξ) dξ and z -(x) = Z R µ -(x, ξ)ψ -(x, ξ) dξ, (3)
where both ψ + and ψ -store a part of the history of the input data u. They are respectively solution to the forward and backward ordinary differential equation in x,

∂xψ + (x, ξ) + θ + (ξ)ψ + (x, ξ) = u(x) with ψ + (0, ξ) = 0, (4)
and

∂xψ -(x, ξ) -θ -(ξ)ψ -(x, ξ) = u(x) with ψ -(1, ξ) = 0, ( 5
)
where ξ is a real parameter. We emphasis that they are independently of a specific operator P . Conversely, the coefficients µ + and µ -, called diffusive symbols, depend on P but not on u. The functions ξ → θ + (ξ) and θ -(ξ) parameterize two closed paths in the complex plane, satisfying the cone condition, and enlacing the singularities of the Laplace transform P ± defined hereafter. The diffusive symbol derivation requires several steps. The two functions y → p(x, x -y) and y → p(x, x + y),

corresponding to the causal part and the anti-causal parts of the impulse response, are analytically extended to R + . We assume that their Laplace transforms P + and P -are well-defined in C + , and that they admit holomorphic extensions vanishing at infinity. The Diffusive symbols are then given by

µ ± (x, ξ) = - θ ± (ξ) 2iπ P ± `x, θ ± (ξ) ´. (7)

III. APPROXIMATIONS AND ALGORITHMS

We introduce the variational formulations for p + and p -and we approximate their solutions by using convenient Galerkin base. The paths -θ ± are chosen so that to enlace, by the right, the poles of Laplace transforms of the basis functions, and the approximations µ N ± of µ ± are then given by the counterpart of the formula [START_REF] Weideman | Parabolic and hyperbolic contours for computing the Bromwich integral[END_REF] applied to the Laplace transform of the impulse response approximations. Once the diffusive symbols are available, a numerical approximation of the integrals in (3) must be derived. The diffusive realization z(x) = z + (x) + z -(x) is evaluated at some points (xn)n=0,..,N ∈ (0, 1), thanks to a trapezoidal rule by

zn = h ξ M X k=-M µ N + n,k ψ + n,k + µ N - n,k ψ - n,k , (8)
for a M ∈ N * , where h ξ > 0 is the step size (distance between two consecutive integration point), µ N n,k and ψ n,k are some discrete values of µ N and ψ at the points xn and ξ k . Such quadrature rule is similar to this encountered in the inverse Laplace transform in the paper of J. A. C. Weideman and L. N. Trefethen [START_REF] Weideman | Parabolic and hyperbolic contours for computing the Bromwich integral[END_REF]. These authors have optimized the parameters of two different classes of contours, typically a parabola

-θ ± (ξ) = θp (iξ + 1) 2 for ξ ∈ R, (9)
and a hyperbola

-θ ± (ξ) = θ h (1 + sin (iξ -α)) for ξ ∈ R, (10)
for some positive real numbers θp, θ h and α the hyperbola asymptotic angle. Our computation are based on the same optimized paths.

In the FPGA, the data are allocated to processing units according to a mapping from the set of space locations {xn | n = 0, .., N } to the set of processing units. For real time computation, the diffusive symbols are computed in a preprocessing step, and only the computation of ψ + n,k and of ψ - n,k , and then of zn through (8), remain to be done at each time step. The summation ( 8) is clearly a local (in space) operation, and the resolution of the Cauchy problems [START_REF] Yakoubi | Deux Méthodes d'Approximation pour le Contrôle Optimal Semi-Décentralisé des Systèmes Distribués[END_REF][START_REF] Hélie | Criterion design for optimizing low-cost approximations of infinite-dimensional systems: towards efficient real-time simulation[END_REF] requires operations between neighbors nodes only. Indeed, we use the forward and backward recurrence relations associated to a discretization method with piecewise constant interpolation of u,

ψ + n+1,k = ψ + n,k e -θ + (ξ k )(x n+1 -xn) + e -θ + (ξ k )(x n+1 -xn) -1 -θ + (ξ k ) u(xn), with ψ + 0,k = 0, and ψ - n,k = ψ - n+1,k e -θ -(ξ k )(x n+1 -xn) - e -θ -(ξ k )(x n+1 -xn) -1 -θ -(ξ k ) u(xn), with ψ - 1,k = 0, (11)
easy to establish, see [START_REF] Yakoubi | Deux Méthodes d'Approximation pour le Contrôle Optimal Semi-Décentralisé des Systèmes Distribués[END_REF]. Clearly, the approximation ψ + n+1,k (resp. ψ - n,k ) at a node xn+1 (resp. xn) is derived from the approximation ψ + n,k (resp. ψ - n+1,k ) at the previous node xn (resp. next node xn+1) for all k. Finally, we present the algorithms in a compact form, where

h = xn+1 -xn is the spatial discretization step, α ± k = e -θ ± k h , β ± k = α ± k -1 -θ ± k , γ ± k = α ± k -θ ± k , and θ ± k = θ ± (ξ k ) for convenient discretization points ξ k . Algorithm 1 Diffusive Realization of z + (x) 1: Offline Computation of diffusive symbol µ N + (x, ξ) 2: Online Computation 3: for n = 0, ..., N do 4: for k = -M, ..., M do 5: ψ + n+1,k = α + k ψ + n,k + β + k un, ψ + 0,k = 0, 6:
end for 7:

z + n+1 = h ξ P M k=-M µ N + n+1,k " α + k ψ + n,k + γ + k un " 8: end for Algorithm 2 Diffusive Realization of z -(x) 1: Offline Computation of diffusive symbol µ N -(x, ξ) 2: Online Computation 3: for n = 0, ..., N do 4:
for k = -M, ..., M do 5:

ψ - n,k = α - k ψ - n+1,k -β - k un, ψ - N ,k = 0, 6:
end for 7:

z - n = h ξ P M k=-M µ N - n,k " α - k ψ - n+1,k -γ - k un " 8: end for

IV. SCALING AND VALUE ENCODING

In this section we explain how we proceed to prepare Algorithm 1 before its implementation in a FPGA. Since Algorithm 2 is very similar to that one we only considered the former.

With the MATLAB implementation of Algorithm 1, all values are encoded with floating point type. More precisely, the MATLAB code uses complex numbers encoded with floating points. By default, the FPGA uses only a bit vector for coding integers and the use of floating points requires many logical cells on the FPGA. So, for many kinds of FPGA, the use of floating is not recommended. In this work, in order to be efficient, one of our constraints was to compute only with integers. In order to scale the values of the variables involved in the computation, we started by analyzing the algorithm. Doing this, we observed that the following variables could be scaled: β, γ, µ and u. Since β and γ depend on α, we do not scale alpha which is almost naturally scaled without doing anything. So, for γ, µ and u, we compute the maximum values, respectively called γ max , µ max and umax. Then we compute the scaled variables, respectively called β scaled , γ sclaed , µ scaled and u scaled with β scaled = β/γ max , γ scaled = γ/γ max , µ scaled = µ/µ max and u scaled = u/umax. It should be noticed that as β and γ depend on α we chose to scale those variables by γ max which is larger than β max .

So in Algorithm 1, we can deduce that ψ (line 5) is automatically scaled since its computation involves scaled variables. Likewise, the computation of z (line 7) is scaled and involves the products of variables scaled by µ max , γ max and umax so in order to obtain the final result, we need to multiply the result by umax ×

β max × µ max . Thus, z = z scaled × umax × β max × µ max
To be efficient on FGPA architecture, one must avoid, when possible, floating point computation. So, next we tried to convert floating point values into integer. Using integers encoded on 32 bits, results of the computation are very similar to the ones using floating points. Nevertheless, the smaller the number of bits used to code integer in the FPGA architecture is, the less numerous are the logical cells used, so the more parallelism capabilities we obtain.

That is why, next we made some experiments using MATLAB to determine the influence of the number of bits to code integer. During the computation some variables take values greater than one even if they were all scaled before the computation. So we had to measure the number of bits required to be added in order to take this into account. In practice, for the experiments we performed, only two bits were sufficient. Likewise, as the FPGA is only ableto manipulate vector of bits representing positive integer, we had to add one more bit for negative numbers. In Figure 1, we compare the results for the causal diffusive realization with floating points, with integers encoded on 16 bits and with the exact solution. The three curves seem very similar even with the computation using integers. In Table I In Figure 2, we highlight the most significant errors due to the influence of the small number of bits used to encode integers with the diffusive representation. We can see the observation highlighted in Table I concerning the number of bits.

V. FPGA FPGA unlike to CPU has no predefined functionalities. This allows to have a finer granularity to implement algorithms.

The FPGA used for this work is a Xilinx Spartan3A ( [START_REF]Spartan 3A datasheet[END_REF]) with: Several implementations of a given algorithm are possible thanks to the flexibility of this type of component:

• a classical iterative implementation, as done in a general purpose computer. This solution is not efficient because for many treatments, waiting the last processing step to be completed before starting the next is needed, but the FPGA gate area use is low, • a pipeline solution. With this strategy, the next treatment can be started on the next clock event. This type of implementation exhibits better global time performances since each new treatment only adds one clock period, but requires more gate area than the iterative technic since registers are needed to store temporary results. • parallel processing. Many independants treatments could be made in the same time. This type of implementation has the advantage to reduce the global time and the critical time is equal to the longest treatment.

VI. POSSIBLE SOLUTIONS OF IMPLEMENTATION TO FPGA

In this section, we explain how we implemented 1 on an FPGA. The implementation choices are directly linked to constraints of the equations (variable dependencies) and FPGA capabilities.

Figure 3 illustrates the relationships between ψ, u(xn), and ξ k . As presented in Figure 3, for a same ξ, ψ(xn+1, ξ) is dependent on ψ(xn, ξ) value. It is not possible to start the next calculation step before having finished the previous one. However, the computation of ψ for different ξ are independent.

Figure 4 shows the relation between the contribution of P (un): there are fewer constraints than ψ because it is possible to add a new value to an existing one or make every addition before to store the result.

The other aspect to be considered is the FPGA ressource availability. For example, if an equation needs to do 50 multiplications with a FPGA containing only 16 multiplier blocks, pipelining the computation to be impossible. The pipeline needs to have at least as many multipliers as required multiplications. The same problem arises for memory management (RAM) and multiplier, a scarce resource on FPGA. These constraints are important because they guide choices in terms of possible solutions, and the manner to implement them. For example, because the value of ψ(xn, ξ k ) is needed for start the computation of ψ(xn+1, ξ k ), it is not possible to use a pipeline solution in this direction.

ψ(0, 3) ψ(0, 2) ψ(1, 2) ψ(1, 1) ψ(2, 1) ψ(0, 0) ψ(2, 3) ψ(0, 1) ψ(1, 3) ψ(1, 0) ψ(2, 2) ψ(2, 0) ξ3 ξ1 ξ2 ξ0
In the sequel, we detail the three solutions above introduced. This is not exaust. Only possible solutions according to the constraints of equations and available ressources are presented.

A. Iterative

The first solution consists in iterating on ξ values. For each value of ξ, it consists in computing the ψ recurrence. This solution is economic in logic cells because data accesses are made sequentially, it is possible to use RAM blocks. But due to the relation between every ψ for a specific ξ k , the total computation time is equal to the time needed to perform the computation for one value of ψ, multiplied by the number N of inputs u and the number of lines (iterations needed for the convergence M ).

B. Pipeline

The second solution is to use a pipeline mechanism. To be able to use it, it is necessary to compute every ψ(xn, 0) values, and next to compute ψ(xn, 1), and so on. This solution is better than the first one because the next start of computation is made one clock period after the previous one. As in the first solution, it is possible to use block RAM because data are accessed sequentially. However, this solution has two drawbacks:

• this solution needs to have the same amount of multipliers than multiplications in equations, • it needs Block RAM to store ψ values to make them available for the next ξ.

C. Parallel

The last solution consists in implementing the computation in parallel. Because there is no dependency between ψ(xn, 0) and another ψ for a different ξ, it is possible to achieve every computation on the same time. This solution is better for the computation time because the size of ξ has no effect on the global computation time. However because constants values are needed at the same time, it is not possible to use RAM. Hence, the resource usage is more important than for the two previous solutions, and the size for N and M will be limited by the amount of logical cells available in the FPGA.

In the following, we only consider the last solution.

VII. IMPLEMENTATION

On the one hand, for any couple of values for xn and ξ k , the calculation of ψ requires 15 multiplications. On the other hand, the Spartan3A has only 16 multipliers. Consequently, it is necessary to distribute a subset of multipliers to each computing branch. The consequence is that it is not possible to dedicate a multiplier to each multiplication. Each component needs to be used more than once.

In order to do that, the only solution consists in creating a Finite State Machine (FSM). Equations are split into multiple states, the number of states depending on the amount of multipliers blocks available for a branch. For each state, the output of a multiplier is affected through addition or substraction to registers, and the next computed values are affected to the multiplier inputs.

The complete structure is split into two components:

• the first is dedicated to handle state evolution, steering RAM.

• the second is dedicated to load values on multipliers and to store results.

This structure has been used for following reasons:

• reducing logical cell usage by having only one state machine,

• a RAM or a signal cannot have more than one source. Other solutions would have arisen problems to handle RAMs.

VIII. AREA USAGE OF FGPA

The chosen solution has the drawback of requiring many logical cells to store constants. Consequently, with a Spartan3A FPGA, the maximum size for constants is illustrated in Figure 6 where optimal values seems to be N = 15 and M = 4 or N = 8 and M = 6.

The evolution curve depends on many parameters: In consequence, for any other encoding integer number of bits, the area usage of logical cells is different.

It evolution of area usage as a function of these parameters is difficult to predict, and is here considered experimentally. For example, the increase on area usage is different for 1 bit added for different M or N magnitude. But, increasing N has a lower impact than increasing M , since N changes only the size for u (stored in RAM) and for µ.

On the other hand, increasing M increases the size of every constant stored with logical cells.

IX. RESULTS

In order to validate our FGPA implementation and measure the execution time of each solution, we compare results synthesized on the FGPA, running on MATLAB, and the equivalent program implemented in C and running on a general purpose computer for a parameter set M = 4 and N = 8. Numerical results are exactly identical whether resulting from the execution on the FPGA, on MATLAB or using the C programs: the algorithm implementation is thus considered correct. For comparing elapsed time, we executed the C program on a laptop computer with an x86 processor 1.6GHz. It takes approximately 60, 000ns to run the simulation.

Time results are given in Figure 7: we observe that on the FPGA, the computation time is approximately 720ns, or an improvement of nearly two orders of magnitude with respect to the general purpose computer.

Figure 8 exhibits the extrapolation of the computation time for others solutions. These results were obtained using the following assumptions

• t ψ = nstate * t clock • t parallel = t ψ * M • titerative = t ψ * M * N • t pipeline = t ψ + ((M * N ) -1) * t clock with:
• t ψ = time to compute one occurence of ψ,

• nstate =number of FSM states, • t clock = clock period. For the same values M = 4 and N = 8, we estimate the following computation times:

• iterative solution: 2800 ns, • pipeline solution: 400 ns. for N = 8 and M = 4, the pipeline solution is better. Indeed the advantage of a parallel implementation appears when N < n stat -1 n stat -M with nstat the number of state for solved an occurence (Fig. 9).

X. CONCLUSION

We have shown that the diffusive representation is well suited for FPGA implementation of the realization of a linear operator which plays a central role in a distributed filtering or control law. We have developed every steps to analyze computation accuracy. In particular, we have tested data formats consistent with hardware implementation. Our results proved a high efficiency in terms of computation time. However, using a Spartan3A FPGA, the amount of available values for constants and for inputs u is limited by the number of logical cells. In a future work, we plan to implement an evolution of this control law. In particular, we will focus on the practically useful case where M is lower than N , and implement the pipeline mechanism that would yield best performances.

Fig. 1 .

 1 Fig. 1. Results obtained by the exact solution of the problem and by the two different approximations (one with floating points and the other with integers encoded on 16 bits).

Fig. 3 . 2 Fig. 4 .

 324 Fig. 3. Dependency for ψ computation

4 Fig. 5 .

 45 Fig. 5. Example of FSM to implement equations

Fig. 6 .

 6 Fig. 6. Percentage of logical cells usage with respect to M and N for a 9 bits integer encoding

Fig. 7 .Fig. 8 .Figure 8

 788 Fig. 7. Generation of a signal indicating the begining and the end of the compute, observed using a 60 MHz bandwidth oscilloscope.

Fig. 9 .

 9 Fig. 9. Simulated computation time for M = 20 and 1 ≤ N ≤ 16

  , we can compare the maximum error and the relative error between different integer computations encoded with different numbers of bits and the exact solution. We can see that the best number of bits is 16 and that the second best is 14 instead of what we should have expected with 15.

	Number of bits	Maximum error	Relative error
	16	9.21e-5	1.05%
	15	2.08e-4	2.51%
	14	1.38e-4	1.51%
	13	4.82e-4	6.79%
	12	6.87e-4	9.23%
	11	1.4e-3	18.69%
		TABLE I
	INFLUENCE OF THE NUMBER OF BITS ON THE ERROR COMPARED TO THE
		EXACT SOLUTION.

• the number of significant bit for values, constants and results,• the magnitude of N and M .


Mélanie Jucla

Samuel Planton

Jean-François Démonet

Christiane Soum-Favaro

ERP and behavioral evidence for interaction/cascade between central (linguistic) and peripheral (motor) processes during word handwriting

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

ERP and behavioral evidence for interaction/cascade between central (linguistic) and peripheral (motor) processes during word handwriting.

Mélanie Jucla 1 , Samuel Planton 1,2 , Jean-François Démonet 2,3 & Christiane Soum 1 1 E.A Octogone. Laboratoire Jacques-Lordat, E.A 4156, Université Toulouse II Le Mirail, Toulouse, France 2 Inserm; Imagerie cérébrale et handicaps neurologiques UMR 825; F-31059 Toulouse, France. 3 Centre Leenaards de la Mémoire, CHUV Département des neurosciences cliniques, Lausanne, Switzerland Contact : melanie.jucla@univ-tlse2.fr

Background

Models of handwritten language production make a distinction between central (access to semantic, orthographic and phonological information) and peripheral (allographic and gesture planning) processes (see [START_REF] Van Galen | Handwriting: Issues for a psychomotor theory[END_REF][START_REF] Rapp | The integration of information across lexical and sublexical processes in spelling[END_REF]. Though, this is still a matter of debate whether these central and peripheral modules are processed in a cascaded or in a serial way [START_REF] Delattre | Written spelling to dictation: Sound-to-spelling regularity affects both writing latencies and durations[END_REF][START_REF] Damian | Advance planning of form properties in the written production of single and multiple words[END_REF]. In the same view, another question still debated is "do central and peripheral processes 'interact' in handwritten word production compared to oral naming"?

ERP data Variables measured : Errors, writing latencies and duration, relative writing speed STATISTICS: Linear mixed models (subjects and items as random effects). Length x regularity x position effects RESULTS:

• Regularity effect Latencies ì (Fig. 2) and speed î for words irregular at the beginning • Length effect (controlled for stimulus length) Latencies ì and speed ì for short words Experiment 2

Stimuli : 180, monomorphemic nouns of 1, 2 et 3 syllables. Consistency manipulated on the first and last segment (à 4 conditions, in/consistent at the beginning and/or at the end) RESULTS:

• Consistency effect Latencies ì (Fig. 3) for words inconsistent at the beginning Speed î for words inconsistent at the end

RESULTS and DISCUSSION

We obtained typical 'N400' interference effect in the Stroop task mainly at medial sites in the oral modality. We put forward a different scalp topography in the written modality [START_REF] Perret | Comparison of electrophysiological correlates of writing and speaking: a topographic ERP analysis[END_REF] with an interference effect that only starts around 450 ms. Lateralization is partly due to motor response and thus could partially explain the late interference effect. This effect could account for the specificity of hand-written language production supporting the view that peripheral processing start before lexical and orthographic selection ends and might impact the 'conflict resolution' à In favor of interaction between central and peripheral processes

Figure 1 :

 1 Figure 1 : ERPs acquired in the a) oral and (b) written color Stroop task. Waveforms were obtained from linear derivation of:

•

  Length effect (controlled for stimulus length)Latencies ì for one syllable monomorphemic nouns, varying on length (short -1 syll.-vs. long -3 syll.-) and regularity paired on frequency(Lexique.org). Irregularity at the beginning or at the end

Figure 2 :

 2 Figure 2 : Latencies are longer for words irregular at the beginning

Figure 3 :

 3 Figure 3 : Latencies are longer for words inconsistent at the beginning

Figure 4 :

 4 Figure 4 : Relative writing speed during word writing for long consistent and inconsistent words


Bioindication and Collembola

JF PONGE

Soil springtails are abundant in a variety of environments, from forests to agricultural crops and bogs, from the deep soil to the crown of trees and from seashore to higher mountains [START_REF] Ponge | Biocenoses of Collembola in atlantic temperate grass-woodland ecosystems[END_REF][START_REF] Sadaka | Soil animal communities in holm oak forests: influence of horizon, altitude and year[END_REF][START_REF] Mouloud | Fine scale changes in biodiversity in a soil-water ecotone: Collembola in two peat-bogs of Kabylia (Algeria)[END_REF]. If most species prefer moist environments, due to their tegumentary respiration, some species are adapted to desert and other psammophilic environments, where they live in small interstices between sand grains [START_REF] Greenslade | Survival of Collembola in arid environments: observations in South-Australia and the Sudan[END_REF][START_REF] Thibaud | The Collembola of littoral sands of Madagascar[END_REF], and many collembolan species support saline stress quite easily [START_REF] Witteveen | Osmotic and ionic regulation in marine littoral Collembola[END_REF][START_REF] Owojori | Comparative study of the effects of salinity on life-cycle parameters of four soil-dwelling species (Folsomia candida, Enchytreus doerjesi, Eisenia fetida and Aporrectodea caliginosa)[END_REF]. However, this apparent ubiquity masks profound differences in the ecological requirements of species. As a consequence, springtail communities exhibit a variable composition. The most cited examples are those relative to soil acidity. It has been shown that, beside a basic species pool, communities of acidic (ph < 5) and less acidic to neutral or alkaline environments differ markedly at the species level [START_REF] Hågvar | Collembola in Norwegian coniferous forest soils. III. Relations to soil chemistry[END_REF][START_REF] Ponge | Biocenoses of Collembola in atlantic temperate grass-woodland ecosystems[END_REF][START_REF] Loranger | Does soils acidity explain altitudinal sequences in collembolan communities?[END_REF]. Although the direct sensitivity of Collembola to pH has been experimentally demonstrated [START_REF] Salmon | Ionic identity of pore water influences pH preference in Collembola[END_REF], other studies showed that the attractiveness of some springtails to earthworms and the refuges they find from predation in networks of earthworm galleries may explain why some species are absent from acidic soils deprived of burrowing animals [START_REF] Salmon | Earthworm excreta (mucus and urine) affect the distribution of springtails in forest soils[END_REF][START_REF] Salmon | Earthworms and Collembola relationships: effects of predatory centipedes and humus forms[END_REF]. Agricultural soils differ also from forest environments, probably due to a complex and still unexplained factor embracing both climate and soil conditions [START_REF] Ponge | Collembolan communities as bioindicators of land use intensification[END_REF][START_REF] Auclerc | Exprimental assessment of habitat preference and dispersal ability of soil springtails[END_REF]. It has been shown that the first stadium, sometimes called the larval stage, may exhibit specific requirements, not allowing some species to tolerate the contrasted climate of open environments [START_REF] Betsch | Caractérisation de deux phases juvéniles d'Allacma fusca (Collembola, Symphypleona) par leur morphologie et leur écophysiologie[END_REF]. In sub-tropical and Mediterranean areas Collembola may totally disappear during drier months [START_REF] Bandyopadhyaya | Effects of some physical factors and agricultural practices on Collembola in a multiple cropping programme in West Bengal (India)[END_REF], passing the unfavourable season in the form of diapausing eggs or other forms of resistance [START_REF] Poinsot-Balaguer | Comportement des microarthropodes du sol en climat méditerranéen français[END_REF]. For all these reasons, springtail communities have been shown to be good indicators of environmental health [START_REF] Van Straalen | Community structure of soil arthropods as a bioindicator of soil health[END_REF], although it can occur that the influence of site is higher than that of vegetational types, as this has been shown to occur in Egypt [START_REF] Alassiuty | The influence of vegetational type on seasonal abundance and species composition of soil fauna at different localities in Egypt[END_REF], pointing to the need for a comprehensive exploration of the regional pool before reaching conclusions based on species composition. The existence of communities locally depauperated in species (compared to the regional pool) may indicate abrupt changes in land use [START_REF] Ponge | Decreased biodiversity in soil springtail communities: the importance of dispersal and landuse history in heterogeneous landscapes[END_REF] or pollution events [START_REF] Gillet | Changes in species assemblages and diets of Collembola along a gradient of metal pollution[END_REF] with which communities were unable to cope. Similarly, soil springtail communities are sensitive to pesticides and other chemicals in use in agriculture, and their paucity in species may indicate negative impacts on non-target organisms [START_REF] Wiktelius | Effects of insecticides on non-target organisms in African agroecosystems: a case for establishing regional testing programmes[END_REF][START_REF] Ponge | Interaction between humus form and herbicide toxicity to Collembola (Hexapoda)[END_REF].


Anne Verroust
email: anne.verroust@inria.fr

Francis Lazarus
email: lazarus@sic.univ-poitiers.fr

M.-S.I I R C O C Sp2mi

Extracting Skeletal Curves from 3D Scattered Data

We introduce a method for extracting skeletal curves from an unorganized collection of scattered data points lying on a surface. These curves may have a tree like structure to capture branching shapes such as blood vessels. The skeletal curves can be used for different applications ranging from surface reconstruction to object recognition.

Introduction.

The description of objects with cylindrical shapes has been extensively used in the computer vision community [START_REF] Agin | Computer description of curved objects[END_REF][START_REF] Soroka | Generalized cylinders from local aggregation of sections[END_REF][START_REF] Shani | Splines as embeddings for generalized cylinders[END_REF][START_REF] Zerroug | Volumetric descriptions from a single intensity image[END_REF]. In a first approximation it is often convenient to model a shape as a generalized cylinder. This is particularly relevant in the context of reconstruction of anatomical shapes from CT data or scanner [START_REF] Pisupati | A central axis algorithm for 3D bronchial tree structure[END_REF][START_REF] Williams | Form from function: a vector field based approach to the analysis of ct images of the vascular tree[END_REF][START_REF] Näf | Characterisation of 3D organ shapes in medical image analysis using skeletonization[END_REF]. Numerous works have been done for object description using tree-like structures [START_REF] Sequeira | Three-dimensional modeling of tree-like anatomical structures[END_REF][START_REF] Nazarian | Interactivity and Delaunay triangulation for the reconstruction of tubular anatomical structures[END_REF][START_REF] Nazarian | Automatic reconstruction of irregular tubular structures using generalized cylinders[END_REF] or cylindrical representations [START_REF] Burdin | Modeling and analysis of 3-D elongated shapes with application to long bone morphometry[END_REF]). What makes the use of cylindrical descriptions and generalized cylinders so popular is their relative simplicity. A generalized cylinder is parameterized by a spin curve together with a set of cross sections [START_REF] Agin | Computer description of curved objects[END_REF]. Depending on the application, the recovery of such parameters usually requires less work than it would be for a more general description of a shape.

Note that a spin curve can be used as a skeleton of an implicit surface model [START_REF] Bloomenthal | Interactive techniques for implicit modeling[END_REF][START_REF] Bloomenthal | Convolution surfaces[END_REF]: the surface of the object is defined as an iso-surface of a field function generated by the skeleton. In [START_REF] Bittar | Automatic reconstruction of unstructured 3D data: Combining a medial axis and implicit surfaces[END_REF][START_REF] Ferley | Skeletal reconstruction of branching shapes[END_REF] the authors propose an implicit reconstruction of shapes from a set of 3D scattered points using the Voronoï graph of the data points to built a geometric skeleton. Our method is an alternate skeleton construction. It avoids the computation of the 3D Voronoï graph.

Other applications start with volume images, so that the traditional medial axis transform can be applied [START_REF] Näf | Characterisation of 3D organ shapes in medical image analysis using skeletonization[END_REF]. As opposed to voxelized data, we start with an unorganized cloud of points sampled on a surface. A 3D curve must be specified inside the cloud of points in order to use a cylindrical model and to represent the surface. When the set of 3D points is given as a set of 2D planar scans, Burdin et al. [START_REF] Burdin | Modeling and analysis of 3-D elongated shapes with application to long bone morphometry[END_REF] use the curve defined by the centroids of each 2D section. However, this implies a consistent orientation of the cross-sectional planes relatively to the position of the object. Figure 1 shows two different orientations of the cross-sectional planes. Clearly, configuration (a) is not adapted. For some objects, because of their bend shape , it is even impossible to find a suitable orientation for the sectional planes.

Although essential, the issue of finding a suitable axis in order to design a generalized cylinder was rarely discussed. In practical cases, this task can be tedious as mentioned by Nazarian & al. [24]. These authors propose two methods to automate the determination of the axis: the first one consists in a recursive subdivision of the set of data points. At each step the subsets of data points are split by a plane perpendicular to their main axis of inertia and passing through their barycenter. The resulting axis is a polyline composed of main axes of parts of the initial set of data points (cf. Figure 2.a). The authors recognize their method does not always work and show a counterexample where the main axis of inertia, represented in dashed line in Figure 2.b, cannot be used to find the central axis.

the second method works with a set of data points given by planar 2D scans. It uses the triangulation computed by Boissonnat's method [START_REF] Boissonnat | Shape reconstruction from planar cross sections[END_REF] as an approximation of the tubular shape.

In this paper, we present a method that automatically computes a set of 3D curves in a cluster of points. When the points lie on a surface these curves may be used as axes for a cylindrical approximation of the surface. These curves may form a tree structure if the surface resemble a "tree of generalized cylinders". The method works independently of the bending of the surface. As a consequence it will return a correct result for the shape shown on Figure 2.b. Like in other reconstruction algorithms [START_REF] Attali | r-regular shape reconstruction from unorganized points[END_REF][START_REF] Mencl | Graph-based surface reconstruction using structures in scattered points sets[END_REF][START_REF] Amenta | A new voronoi-based surface reconstruction algorithm[END_REF], we make some hypothesis on the input cluster of points.

1:

The data points should be sampled on a certain (unknown) surface.

2:

The sampling should be done in such a way that the distance between neighbor points is small with respect to the width of the tubular parts of the surface. The computation of the axes is based on a decomposition of the surface to be reconstructed into level sets of a scalar function. In [START_REF] Lazarus | Courbes, cylindres et métamorphoses pour l'image de synthèse[END_REF][START_REF] Lazarus | Décomposition cylindrique de polyèdre et courbe squelette[END_REF][START_REF] Lazarus | Level set diagrams of polyhedral objets[END_REF] we used a similar decomposition applied to polyhedral surface. In this work a simplified discrete version of the Morse theory was developed (see [START_REF] Shinagawa | Surface coding based on Morse theory[END_REF][START_REF] Kergosien | Topology and visualization: from generic singularities to combinatorial shape modelling[END_REF][START_REF] Takahashi | Algorithms for extracting correct critical points and constructing topological graphs from discrete geographical elevation data[END_REF][START_REF] Hart | Morse theory for implicit surface modeling[END_REF] for use of Morse Theory in Graphics).

In the present case, we are dealing with a cluster of points and no topological information, such as adjacency relationship, is available. We cannot apply the method of [START_REF] Lazarus | Courbes, cylindres et métamorphoses pour l'image de synthèse[END_REF][START_REF] Lazarus | Décomposition cylindrique de polyèdre et courbe squelette[END_REF][START_REF] Lazarus | Level set diagrams of polyhedral objets[END_REF] which uses this information to build the decomposition of the surface.

Here, our goal is to compute the axes rapidly, rather than reconstructing the surface. We follow the main steps of [START_REF] Lazarus | Courbes, cylindres et métamorphoses pour l'image de synthèse[END_REF][START_REF] Lazarus | Décomposition cylindrique de polyèdre et courbe squelette[END_REF] and adapt them to our particular case. The scalar function used in the decomposition is an approximation of a geodesic distance function to a given point. This function is determined by a single point from the set of data points called the source point. The source point may be interactively designated by the user. It may also be automatically selected using a simple heuristic as explained at the end of Section 3. Intuitively, the source point corresponds to a polar extremity of the shape as shown on Figure 3 decomposition into level sets of a scalar function is done, a single axis is obtained for each cylindrical part associated to it. .

Overview

The computation of the 3D curves is made in four steps, as illustrated in Figure 4. Given a set of data points,

1:

We build the neighborhood graph of order n. This is a non oriented graph whose vertices are the input data points.

The edges have the form (v,v 0 ), where the vertex v 0 or the vertex v is one of the n nearest neighbors of respectively v or v 0 . The subgraphs of the neighborhood graph induced by a point and its neighbors may be non planar and, in general, it will not be possible to fit a suitable polygonal surface to this graph.

2:

From this neighborhood graph and a source point, we construct a distance map and an oriented subgraph called the geodesic graph. The geodesic graph is obtained by computing shortest paths of the neighborhood graph between all vertices and the source point. The distance map is simply the length of these shortest paths. The construction of the neighborhood graph, the distance map and the geodesic graph are described in details in Section 3.

3:

We further compute k level sets of the distance map using the neighborhood graph, where k is a user defined value. The level sets are composed of points at a constant distance from the source. We connect these points on a proximity criterion and identify the connected components. We also compute the centroids of each of these connected components.

4:

The axes are finally obtained by connecting the centroids of successive levels as previously computed. When the levels contain several connected components, the connection of the centroids are based upon the paths of the geodesic graph.

The Geodesic Graph

We are given a set of data points X = fx 1 ; ::::; x m g sampled on an unknown surface S and an integer n. The neighborhood graph N (X;n) is an undirected graph defined as follow:

1: the vertices of N (X;n) are the points of X 2: a pair (x i , x j ) is an edge of N (X;n) if either x i is one of the n nearest points of x j in X or x j is one of the n nearest points of x i in X.

In order to compute the n nearest points of each vertex, we use a subdivision of the space into voxels. We subdivide the enclosing box of the data points into a regular grid of voxels and associate each data point with the voxel containing it. The search of the n nearest neighbors of a data point x is made in two steps. We first visit the voxels close to the one containing x progressively moving away until n points (v 1 ,...v n ) are found. Then, we put d = max 1 i n (distance (v i ; x)) the maximum distance between x and these n points. As we have already found n points (v 1 ,...v n ) located at a distance lower or equal to d from x, we are sure that the n nearest neighbors of x belong to the sphere S(x; d) with center x and radius d. We search all the data points belonging to S(x; d) and we extract the n nearest points on this set, using a quick sort on the distance values.

The number n should be chosen so that N (X;n) is a connected graph. This choice depends on the regularity of the sampling step of the data points on the surface. In the example presented in Figure 8 we used n = 5. When the data points are not uniformly spread over the surface, larger values may be necessary. For instance, in the case of 2D planar scans, where the resolution is smaller inside the scan planes than between the cross-sections, we found that a value of n = 10 was necessary to obtain a connected graph. This was the case for the blood vessels shown in Figure 10.

Together with the selected source point x s , the neighborhood graph is used to compute the geodesic graph G (X;n;xs) .

The geodesic graph is composed of geodesic paths of the neighborhood graph emanating from the data points and joining the source point. The length of an edge of the neighborhood graph is taken as the Euclidean distance between its endpoints. We use Dijkstra's algorithm [START_REF] Cormen | Introduction to Algorithms[END_REF] to efficiently compute these paths and their length. Precisely, Dijkstra's algorithm compute the geodesic distance from any point to the source as well as a pointer to the next point in a geodesic path joining the source. Note that there may exist several geodesic paths from a point to the source. We simply use the path returned by Dijkstra's algorithm. The distance map introduced in Section 2 is defined for every data point to be its geodesic distance to the source. It is extended to the interior of the edges of the geodesic graph with a linear interpolation of the endpoint distance values of the edges.

As we supposed that the edge of the neighborhood graph should lie on the surface, the geodesic paths previously computed are intended to approximate geodesic paths on the surface. Note that The computation of geodesic paths for polyhedral surfaces have been the subject of numerous publications [START_REF] Mitchell | The discrete geodesic problem[END_REF][START_REF] Opitz | Computing shortest paths on polyhedra: applications in geometric modeling and scientific visualization[END_REF]. However, we cannot apply those methods since our data points do not have a polyhedral structure.

Note that G (X;n;xs) is a connected subgraph of N (X;n) since all the geodesic paths share the source point. G (X;n;xs) is actually a tree since Dijkstra's algorithm selects only one path from any vertex to the source. An example of a geodesic graph is shown Figure 6 for a set of data points registered on a femur1 .

The source point x s may be selected by the user or automatically computed by the following two step heuristic:

1 a point x 0 is chosen at random (e.g. the first point in the data set).

2 x s is defined as the furthest point from x 0 on the neighborhood graph.

When the object has a "tubular structure", even with branches, the furthest point from x 0 is a point located at the extremity of a tubular part of the surface. This point can be obtained as the last visited vertex when running Dijkstra's algorithm a first time with source point x 0 .

Levels and branches

After computing the geodesic graph and the distance map, we extract k level sets. Intuitively, each level set corresponds to a set of "sections" of the decomposition of the surface given by the set of data points X = fx 1 ; ::::; x m g.

The number k is given by the user to specify the number of points used to approximate the longest axis. In practice k can be chosen as the ratio of the distance of the furthest point to the average edge length. The k levels d 1 ,. . . , d k are uniformly distributed over the interval d max ; d max ],

where d max is the extremum of the distance map and and are two percentage coefficients. We used = 3% and = 97% on Figure 11. The ith level set S i is the set of points on the neighborhood graph whose (extended to edge interiors) distance is d i .

Next, we need to partition each level set S i into subsets corresponding to the different branches of the surface.

Here, we cannot proceed as in [START_REF] Lazarus | Courbes, cylindres et métamorphoses pour l'image de synthèse[END_REF][START_REF] Lazarus | Décomposition cylindrique de polyèdre et courbe squelette[END_REF] where each vertex of the polyhedron can be classified as regular or critical according to the configuration of distances of its neighbors (see Figure 5). In the polyhedral case, the decomposition is obtained by cutting through the level sets associated to the critical points. Also, the level sets of the geodesic distance are efficiently extracted using the adjacency between faces of the polyhedron.

In our case this adjacency relationship is not available. As a result, points cannot be classified as regular or critical which would in turn provide the branching points and the partition of each level set corresponding to the different branches. To obtain this partition, we first construct the neighborhood graph N (Si;2) , linking every point of S i with its two closest neighbors. Then, we remove from N (Si;2) the edges longer than a certain length w i . w i is determined using the hypothesis (2) on the data points: the distance between points is small with respect to the width of the tubular parts. In practice w i is computed using the median of the distribution of the distance values of the second near-est neighbor of each point inside S i . Finally, we partition this graph into its connected components. We also compute the centroid of each connected component by taking the barycenter of its points. The skeletal curves are constructed by connecting the centroids of successive level sets. When the level sets have several connected components we use the geodesic paths to connect the centroids: if a level set S i is partitioned into l i connected components (S i;1 ; ::S i;li ), for each j l i , we select a point of S i;j and follow its geodesic path to the source point until an edge containing a point of S i 1 is encountered. This point belongs to the connected component of S i 1 which is either associated to the branch containing S i;j or located at the birth of this branch. Eventually, we may approximate the axes with spline curves as shown for the femur on Figure 6 and7 . We use a simple fitting procedure from [START_REF] Schneider | An algorithm for automatically fitting digitized curves[END_REF][START_REF] Eilers | Smoothing and interpolation with finite differences[END_REF] to compute these splines from the connected centroids.

Results

We applied our skeletonization method on a set of examples presented from Figure 8 to Figure 11:

Figure 8 shows two models of a man and a horse. These models are public domain models found on the Web pages of Cyberware2 . The source point of the man is located at the top of his head. The horse was scanned from many orientations resulting in 210 linear scans which consisted of a series of scans of different regions (views) of the object. The source point is located at the extremity of its jaw.

Here, the branches are correctly located and a tree-like axial structure is built.

The aorta presented in Figure 9 is described by a set of 2D scans. This data comes from the National Library of Medicine Visible Human Male3 . Our method is little sensitive to the location of the source points: in Figure 9.A, two neighbor points and in Figure 9.B, three points located at different extremities of the tubular parts of the surface have been used to compute the axial structures. The difference between their corresponding axial structures is negligible.

We obtained good results with our method for the two blood vessels of Figure 10 4 . These data sets are especially challenging since they contain many vessels close to each other.

The axis of the colon presented in Figure 11 cannot be found using Nazarian's subdivision method [START_REF] Nazarian | Interactivity and Delaunay triangulation for the reconstruction of tubular anatomical structures[END_REF]. It is de-scribed by a set of 2D scans and comes from the National Library of Medicine Visible Human Male. Table 1 shows the CPU time required to compute the axial structure for the models shown Figures 6 to 11. We run non optimized code on a Indy The computing time depends mainly on the number of data points and on the order of neighborhood graph. In fact, the computation of the neighborhood graph is one of the most expensive part of our method.

Conclusion

We have presented a new technique to associate an axial structure to a set of scattered data points. The interaction is reduced to the selection of a source point and can be automated using a simple heuristic. The main steps for the construction of the axes are:

1 the computation of the neighborhood and geodesic graphs, using Dijkstra's algorithm, 2 the computation of geodesic levels and the location of branching nodes. Our method is rapid and robust with respect to the location of the source point. The "quality" of the resulting skeletal curves mainly depends on the appropriate choice of the order of the neighborhood graph.

Our skeletonization algorithm relies on the location of the source point. This means that the axial structure found by the algorithm is not unique for a given object. The axial structure is actually defined for a pair (object, source point).

The question that arises from this observation is how can we extract interesting axis for objects with tubular shape? We should note that the precise definition of a tubular shape is not clear (except for the trivial case of a curve extrusion along an axis). As a consequence there is no precise definition for the axis of a tubular shape. The medial axis, for instance, does not always produce what we would expect from a tubular object. As opposed to our skeletons the medial axis is in general composed of pieces of curves as well as pieces of surfaces with many ramifications [START_REF] Blum | Biological shape and visual science (part I)[END_REF][START_REF] Goldak | Constructing discrete medial axis of 3D objects[END_REF] It can be simplified to obtain a wireframe skeleton [START_REF] Attali | Squelettes et graphes de Voronoi 2D et 3D[END_REF][START_REF] Attali | Computing and simplifying 2D and 3D semicontinuous skeletons of 2D and 3D shapes[END_REF] but this simplification process is not entirely automatic, as some thresholds must be fixed. The whole process may also be expensive as the simplification process works on the medial axis of the object represented by the set of sampled data points.

There are many objects, however, for which the definition of an axial structure is intuitively so obvious that we could draw this axis without any hesitation. Blood vessels are good examples of such objects. Our initial question is now turned into the more pragmatic question: how does this skeletonization algorithm provide an efficient way for computing axis in blood vessels structure? The only point we have to worry about is the selection of the source point since the computation is entirely automatic once a source point is selected. On one hand it is seen in Figure 9.A that the algorithm is little sensitive to the precise location of the source point. This is due to the continuity of the distance function with respect to the source location. On the other hand we have seen Section 3 that the source point can be selected using a simple heuristic. Figure 9.B shows that this heuristic works well in practice: extremities on different branches give rise to nearly identical axis.

Work in progress include attempts to model more precisely the location of the branching nodes. A more robust construction of the neighborhood graph using an adaptive number of neighbors based on the local density of data points is also under consideration.

Figure 1 .

 1 Figure 1. Two sets of 2D sections sampled on a surface.As opposed to voxelized data, we start with an unorganized cloud of points sampled on a surface. A 3D curve must be specified inside the cloud of points in order to use a cylindrical model and to represent the surface. When the set of 3D points is given as a set of 2D planar scans, Burdin et al.[START_REF] Burdin | Modeling and analysis of 3-D elongated shapes with application to long bone morphometry[END_REF] use the curve defined by the centroids of each 2D section. However, this implies a consistent orientation of the cross-sectional planes relatively to the position of the object. Figure1shows two different orientations of the cross-sectional planes. Clearly, configuration (a) is not adapted. For some objects, because of their bend shape , it is even impossible to find a suitable orientation for the sectional planes.Although essential, the issue of finding a suitable axis in order to design a generalized cylinder was rarely discussed. In practical cases, this task can be tedious as mentioned byNazarian & al. [24]. These authors propose two methods to automate the determination of the axis: the first one consists in a recursive subdivision of the set of data points. At each step the subsets of data points are split by a plane perpendicular to their main axis of inertia and passing through their barycenter. The resulting axis is a polyline composed of main axes of parts of the initial set of data points (cf. Figure2.a). The authors recognize their method does not always work and show a counterexample where the main axis of inertia, represented in dashed line in Figure2.b, cannot be used to find the central axis.the second method works with a set of data points given

Figure 2 .

 2 Figure 2. a: recursive subdivision using main axes of inertia (Nazarian & al); the first main axis is shown in dashed line; The main axes of inertia and the perpendicular splitting planes for the second level of recursion are shown in plain lines. b: The first axis of inertia and the perpendicular splitting plane are not suitable for the subdivision algorithm.

Figure 3 .

 3 Figure 3. The source point and the corresponding decomposition.

Steps 3 Figure 4 .

 34 Figure 4. The main steps of our method.

Figure 5 .

 5 Figure 5. the leftmost vertex is regular and the other are critical (the + (respectively ) associated to a vertex indicates that the geodesic distance between the vertex and the source is greater (respectively lower) than the geodesic distance between the source and x).

Figure 6 .

 6 Figure 6. Up: 15 level sets computed on the geodesic graph of a femur. Down: the corresponding axis.

Figure 7 .

 7 Figure 7. Different number of levels used to compute the axis of the femur : 30 (grey curve) and 90 (black curve).

Figure 8 .

 8 Figure 8. 15 level sets and the corresponding skeletal curves on a man and on a horse.

Figure 10 .

 10 Figure 10. Two ramifications of blood vessels.

Figure 11 .

 11 Figure 11. Up: a colon. Down: its axis (two different views).

Table 1 . Statistics for Figures 6 to 11.

 1

	reference	number of number of	order of the	CPU time in seconds
	of figures	data points	levels	neighborhood graph	on an Indy
	6	16781	15	5	27.83
	7	16781	30	5	33.20
	7	16781	90	5	44.83
	8	21500	15	5	52.21
	8 (horse)	48486	15	6	132.61
	9 (A and B)	7688	50	5	40.01
	9 (A)	7688	50	5	39.54
	9 (B middle)	7688	50	5	35.09
	9 (B right)	7688	50	5	35.91
	10 (left)	11807	60	10	28.85
	10 (right)	13582	50	8	45.49
	11	55573	150	5	235.56

data, courtesy of Jai Menon, IBM T. J. Watson Research Center.

http://www.cyberware.com/models/index.html

http://www.mayo.edu/bir/Models/Visible Human Male/VRML/index.html

Data, courtesy of Jacques Feldmar.