Modesto T Lopez-Lopez
email: modesto@ugr.es

Laura Rodriguez-Arco

Andrey Zubarev

Pavel Kuzhir

Larisa Iskakova

Fernando Gonzalez-Caballero

N-like rheograms of concentrated suspensions of magnetic particles

We investigate the rheograms of concentrated suspensions of magnetic particles obtained under imposed shear rate in parallel plate geometry. We show that under magnetic field application the usual trend of the rheogram, i.e., increasing shear stress for the whole range of shear rates, is altered by the appearance of a region in which the shear stress decreases as the shear rate is increased. The existence of this region gives to the rheograms an N-like shape. The two initial regions (pre-yield regime) of these Nlike rheograms present unstable flow, characterized by the oscillation of the shear stress with time for each imposed value of shear rate. We also show that rheograms obtained at different sample thicknesses approximately overlap in the developed flow regime, whereas there is a tendency of the shear stress to increase as the thickness is decreased in the pre-yield regime. This tendency is likely due to the strengthening of preexisting particle structures by compression as the gap thickness is decreased. Finally, we analyze the effect of the applied magnetic field strength, H, and demonstrate that the rheograms

I. INTRODUCTION

Magnetorheological (MR) fluids are suspensions of micron-sized particles of magnetizable materials. These suspensions are characterized by the tunability of their rheological properties by means of magnetic field application, a phenomenon known as the MR effect [START_REF] Bossis | Ferrofluids[END_REF]. Changes provoked in MR fluids by application of magnetic fields are so intense that their mechanical behavior usually shifts from this of a Newtonian liquid in their absence, to this of a plastic material in their presence [START_REF] Bossis | Ferrofluids[END_REF]. Plastic materials are characterized by a minimum value of stress that must be applied in order to initiate their flow, known as the yield stress. In the case of MR fluids the magnitudes of both the yield stress and the viscosity increase when the intensity of the magnetic field is increased. This smart property is the base of many applications of MR fluids in engineering and biomedicine [START_REF] Park | Magnetorheology: materials and application[END_REF].

Rheology of MR fluids has been an active field of research since the eighties of the twentieth century. Main attention has been paid to the change of the characteristic rheological properties (yield stress, viscosity and viscoelastic moduli) with the intensity of the applied magnetic field; i.e., to the MR effect. The aim has usually been the experimental characterization of new MR fluids, and the understanding from the viewpoint of the particle interactions of the physical origins and strength of the observed MR effect. Recently, MR fluids constituted by fiber-like particles were prepared [START_REF] Bell | Magnetorheology of submicron diameter iron microwires dispersed in silicon oil[END_REF][START_REF] López-López | New magnetorheological fluids based on magnetic fibers[END_REF] and their rheology was extensively studied afterwards [START_REF] López-López | Magnetorheology of fiber suspensions. I. Experimental[END_REF][START_REF] De Vicente | Dynamic rheology of sphereand rod-based magnetorheological fluids[END_REF][START_REF] Bell | Impact of nanowires on the properties of magnetorheological fluids and elastomer composites[END_REF][START_REF] Zubarev | N-Like rheograms of suspensions of magnetic nanofibers[END_REF]. These studies revealed some unusual rheological behavior, which had not been shown before by any other kind of MR fluids. This is the case, for example, of a region of negative differential viscosity, characterized by a decrease of the shear stress as the shear rate is increased [START_REF] Zubarev | N-Like rheograms of suspensions of magnetic nanofibers[END_REF]. This kind of phenomena has stimulated the search of other unconventional rheological behaviors in MR fluids.

Very recently, we have found that near the limit of maximum-packing fraction, conventional (based on spherical particles) MR fluids also display a region of negative differential viscosity [START_REF] López-López | Stick-slip instabilities in the shear flow of magnetorheological suspensions[END_REF]. It is well documented that a stable flow is impossible within the shear rate range corresponding to the negative differential viscosity [START_REF] Quemada | Unstable flows of concentrated suspensions[END_REF] and, in fact, we also found that for a given value of the shear rate pertaining to this region, the shear stress oscillated with time, in a way that resembled dry stick-slip motion [START_REF] López-López | Stick-slip instabilities in the shear flow of magnetorheological suspensions[END_REF]. In the present work, we focus on the analysis of the complete flow curves of concentrated MR fluids. As we show below, the region of negative differential viscosity gives to the rheograms an N-like shape, not observed before for other types of MR fluids constituted by spherical particles. We also analyze in the present work the effect of the measuring system gap thickness on the N-like rheograms. Interestingly, there is a tendency of the shear stress to increase as the gap thickness is decreased in the pre-yield regime.

Finally, we analyze the effect of the strength of the applied magnetic field, H, and obtain an increase of the shear stress with H 1.5 for each given value of the shear rate, for the whole region of shear rate values.

II. MATERIALS AND METHODS

For the preparation of the MR fluids we used carbonyl iron powder (Fe-CC, BASF, Germany) as solid phase and mineral oil (Sigma Aldrich, Germany) as carrier liquid. According to the manufacturer, Fe-CC particles have a silica coating and are spherical in shape with a median diameter of 5 µm. We prepared MR fluids with 50 vol.% approximate concentration, following the protocol described in [START_REF] López-López | Yield stress in magnetorheological suspensions near the limit of maximum-packing fraction[END_REF]. For confirmation of some hypothesis that we will discuss in Sec. III, we also prepared MR fluids containing 10 vol.% of Fe-CC particles, by dispersing proper amounts of the powder in mineral oil.

We characterized the rheological properties of our suspensions at 25 ºC. In order to show that the main findings of our work are not the result of a particular experimental setup, we used two different rheometers and three different measuring geometries, as described in what follows.

(i)

Haake Mars III provided with parallel plates with smooth surfaces. If not specified otherwise, data shown below were obtained with a Haake Mars III rheometer and measuring system geometry consisting of a set of parallel plates, 35 mm in diameter, with smooth surfaces. In particular, we used this setup for the study of the effect of gap thickness, which we varied between 0.1 and 1 mm. For each gap thickness, we calibrated the rheometer by using two viscosity oil standards (Brookfield, USA) of 1.930 Pa•s and 6.470 Pa•s, to ensure that the measurements obtained at the different gap thicknesses could be compared safely.

(ii) Haake Mars III provided with parallel plates with magnetic surface. In order to prove that the findings reported here do not qualitatively result as a consequence of the lack of transmission of shear forces to the measuring geometry due to the formation of a slip layer close to the smooth surface of the moving plate, we performed experiments with setup described in case (i) above, by previously coating the working surfaces by a thin magnetic layernote that this approach completely mitigates the problem of slip on the surface of the rheometer plates (wall slip) [START_REF] Laun | Wall material and roughness effects on transmittable shear stresses of magnetorheological fluids in plate-plate magnetorheometry[END_REF]. For this aim we glued a layer of Fe-CC particles to each of the working surfaces of the set of parallel plates. We checked at the end of the rheological measurements that the layers of Fe-CC particles were still glued to the working surfaces. We worked in this case with a gap thickness of 1 mm.

(iii) Physica-Anton Paar MCR 300 provided with parallel plates with serrated surfaces. A general approach to avoid wall slip is the use of measuring systems with milled, serrated, or otherwise roughened or textured surfaces [START_REF] Buscall | The rheology of concentrated dispersions of weakly attracting colloidal particles with and without wall slip[END_REF][START_REF] Gómez-Ramírez | Wall slip phenomena in concentrated ionic liquid-based magnetorheological fluids[END_REF]. We also tested this approach and for this aim we used a rheometer Physica-Anton Paar MCR 300 provided with a measuring system consisting of a set of parallel plates, 20 mm in diameter, with serrated surfaces (see Figure 1). We worked in this case with a gap thickness of 0.35 mm.

For the application of the magnetic field we used solenoids placed co-axially with the measuring geometry, as described in [START_REF] Kuzhir | Shear and squeeze rheometry of suspensions of magnetic polymerized chains[END_REF]. The experimental protocol used for the rheological measurements is reported in details elsewhere [START_REF] López-López | Stick-slip instabilities in the shear flow of magnetorheological suspensions[END_REF]. Briefly, we proceeded as it follows. First of all, we subjected the MR fluid to a linear shear rate ramp in the range 0.1-100 s -1 , for a total duration of 60 s. Then, we applied a magnetic field of a desired strength and we left the MR fluid at rest for 30 s. Finally, we imposed a shear rate (within the range 0.0025-100 s -1 ) and monitored the shear stress response for 30 minutes. When the measurement at a given shear rate was over, we replaced the sample by a fresh one and repeated the whole process for another shear rate.

III. EXPERIMENTAL RESULTS

A. Time evolution of the shear stress

Under magnetic field application, the time evolution of the shear stress of the MR fluid containing 50 vol.% of Fe-CC revealed that a steady value of the shear stress was not reached even at times as high as 30 minutes, for values of the imposed shear rate below a critical threshold (Figure 2a). As observed, after an initial sharp increase of the shear stress, some saw-tooth-like oscillations appeared at low enough values of the imposed shear rate. In a previous work we studied the origin of these oscillations and showed that from a microscopic viewpoint they were the consequence of the periodic failure and healing of the field-induced particle structures [START_REF] López-López | Stick-slip instabilities in the shear flow of magnetorheological suspensions[END_REF]. In more details, we supposed the appearance of a periodic array of gap-spanning cylindrical columns. These columns were tilted by the motion of the upper rheometer plate, and broken at their middle point at some critical tilt angle (corresponding to the maximum of the stress vs.

strain curve). Afterwards, the two halves of the broken aggregates were supposed to be reformed with other broken columns at the first encounter. The temporal period of the stress oscillations was therefore directly related to the space period of the column-like structure. For more details on the microscopic origin of the stress oscillations see [START_REF] López-López | Stick-slip instabilities in the shear flow of magnetorheological suspensions[END_REF].

As discussed in the section materials and methods we checked that the stress oscillations were not the result of a particular instrumental setup. As observed (see inset of Figure 2a), stress oscillations were also obtained for measuring plates with a magnetic layer coating. In addition, we previously reported similar stress oscillations for measurements at imposed shear rate using rheometer Physica-Anton Paar MCR 300 provided with parallel plates with serrated surfaces [START_REF] López-López | Stick-slip instabilities in the shear flow of magnetorheological suspensions[END_REF], which confirms the independence on the experimental setup of our results from the qualitative point of view. Furthermore, note that similar stress oscillations under simple shear for MR fluids consisting of 30 vol.% of iron in silicone oil were very recently reported by other authors [START_REF] Jiang | Stick-slip behavior of magnetorheological fluids in simple linear shearing mode[END_REF]. Surface roughness of the measuring plates was approx. 0.02 microns, whereas iron particles were 3 microns in diameter in this previous work [START_REF] Jiang | Stick-slip behavior of magnetorheological fluids in simple linear shearing mode[END_REF].

Contrarily to concentrated suspensions, MR fluids containing 10 vol.% of Fe-CC particles did not show stress oscillations, but a steady state regime with an approximately constant value of the shear stress achieved at large enough times (Figure 2b). This fact confirms our previous hypothesis [START_REF] López-López | Stick-slip instabilities in the shear flow of magnetorheological suspensions[END_REF] that stress oscillations observed for highly concentrated MR fluids were a consequence of the stick-slip motion, provoked by the periodic fracture and reformation of the thick 3-dimensional structure induced by the applied field. Note that, for some reason, collective breaking and reformation of aggregates did not occur in suspensions of moderately low particle volume fraction (e.g., 10 vol.%). This was likely because of irregular spacing between columns, which would hinder the collective behavior more easily at low concentrations (when the distance between aggregates was large) than at high concentrations.

Note at this point, that during the decreasing part of the stress oscillations, the imposed value of the shear rate was not maintained (Figure 3). As observed, during the increasing part, the shear rate was perfectly constant and the shear stress monotonically increased. This increase of the shear stress corresponded with the deformation of the 3dimensional particle structure induced by the field. Deformation of this structure required a work, made by the shearing force. At the maximum of the stress within each oscillation, the 3-dimensional particle structure fractured -experimental evidence, as well as detailed discussion, of the microscopic origin of the stress oscillations can be found in our previous work [START_REF] López-López | Stick-slip instabilities in the shear flow of magnetorheological suspensions[END_REF]. Slip of one part of the fractured structure over the other part required much lower shear stress than deformation of the field-induced particle structure and, consequently, the shear rate increased sharply during a short time interval due to the inertia of the upper plate of the rheometer. We refer here to slip occurring between broken aggregates, not to slip between aggregates and the rheometer wall (wall slip). Note that the rheometer was unable to react immediately to the new situation (slip between broken structures) that required much lower shear stress and, consequently, for a short time interval the stress applied provoked a sharp increase of the rotation speed of the upper plate (and thus of the shear rate), as observed in Figure 3. The augmented values of the shear rate elapsed the total duration of the decreasing part of the stress oscillation, presumably until the particle structures reformed again, which marked the onset of the next increasing shear stress region. Both the peak value of the shear rate (with respect to the imposed value), and the duration of the decreasing part of the stress oscillation, decreased with the increase of the imposed value of the shear rate (Figure 3). For the values of applied magnetic field used in our work, the shear rate,  , was maintained within a maximum deviation of approx. 10 -20 % of its imposed value for 0.5

 

s -1 for rheometer Haake Mars III and for 0.025

 

s -1 for rheometer Physica-Anton Paar MCR300 -obviously, this error was much smaller the higher the value of the imposed shear rate, being negligible close and above the yield point.

Such sharp increase of the shear rate was not observed for moderately diluted MR fluids. For these suspensions the shear rate was approximately maintained around the imposed value, as exemplified in Figure 4.

B. N-like shaped rheograms

The analysis of the segments of the stress vs. time curves (like those shown in

Figure 2) at large elapsed times and for each value of the imposed shear rate allowed us to obtain the steady-state rheograms (curves of shear stress vs. shear rate) of our suspensions. For those experimental conditions that gave rise to oscillating values of the shear stress, we took the mean value of the stress (within its oscillations) as the representative stress for the corresponding shear rate -note at this point that if restricted to the shear stress values corresponding to the increasing part of the oscillations, for which the shear rate was perfectly maintained, the same rheograms were obtained.

Concerning the effect of measuring system gap on the rheograms of concentrated MR fluids (50 vol.%), in the absence of applied magnetic field rheograms obtained at different gap thicknesses approximately superimposed for the whole range of shear rates (Figure 5a). On the contrary, when a magnetic field was applied, there were differences in the rheograms obtained at different gap thicknesses (Figure 5b). As observed, in the pre-yield regime (initial increasing part of the shear stress and decreasing part of the shear stress) there was a tendency of the shear stress to increase as the gap thickness was decreased. On the other hand, in the flow regime (shear rate higher than 1 s -1 ) all the curves approximately superimposed.

The enhancement of magnetic field-induced yield stress as the gap thickness was decreased was previously observed for concentrated MR suspensions subjected to compression [START_REF] Tang | Structure-enchanced yield stress of magnetorheological fluids[END_REF][START_REF] Tao | Super-strong magnetorheological fluids[END_REF]. In those cases, compression of pre-existing structures gave rise to thicker particle structures, which resulted in stronger MR effect.

Something similar could be the cause for the enhancement of the shear stress in the pre-yield regime as the gap thickness was decreased, observed in our experimentsnote that in the flow regime particle structures are broken. To check the validity of this hypothesis, we submerged the rotating plate of the Haake Mars III rheometer in a bath (cup) containing the concentrated MR fluid (50 vol.% of iron particles), applied a magnetic field of 5.4 kA/m, and fixed a shear rate of 0.05 s -1 . Under the applied magnetic field and continuous rotation of the plate, we decreased the gap, from 1 mm down to 0.4 mm, stepwise, and monitored the evolution of the shear stress and the normal force (Figure 6).

As observed, the normal force increased abruptly when the gap thickness was decreased, although it relaxed afterwards (Figure 6). The shear stress also increased abruptly when the gap thickness was decreased and then experienced some oscillations in coherence with observations of Figure 2. Note that relaxation of the normal force was very sharp (stepwise) at some specific points that corresponded to the abrupt decrease of the shear stress and, thus, to the breakage of the particle structures into two parts. This result corroborates the hypothesis that under an applied magnetic field, compression results in stronger MR structures. Nevertheless, in the experiments of Figures 2345we compressed our samples without shear and without magnetic field application.

To check if the same phenomenon of structure enhancement could take place under these conditions, we measured the normal force on the upper, rotating plate, submerged in the concentrated MR fluid, contained in the cup, this time without magnetic field and without shear application. Results show that the normal force increased abruptly when the gap thickness was decreased, and relaxed smoothly with time afterwards (Figure 7). The behavior observed in Figure 7 can be explained as it follows. In our concentrated suspensions, very likely there were preexisting particle aggregates as a consequence of colloidal interactions and remnant magnetization. When the gap thickness was decreased, these aggregates would have been compressed, which would explain the observed increase of the normal force (Figure 7). Such compressive pre-load at zero field would have resulted in the interpenetration of asperities of touching particles. Once the field was applied, these particles would have experienced a stronger magnetic attraction than those for which compression was not performed. The lower the gap, the stronger the compression of initial aggregates as inferred from Figure 7 and, therefore, the stronger would have been the interpenetration of the particle asperities, giving rise to a stronger shear stress under magnetic field in the pre-yield regime, in agreement with results of Figure 5b.

Effect of applied magnetic field strength

We also analyzed the effect of the applied magnetic field on the rheograms and, as observed (see Figure 8a), we obtained similar N-like shaped rheograms for different applied field strengths, with a tendency of the shear stress to increase as the magnetic field strength was increased, as expected for a MR fluid. It is worthy to mention at this point that in a previous paper [START_REF] López-López | Yield stress in magnetorheological suspensions near the limit of maximum-packing fraction[END_REF] we developed a model for the magnetic field-induced yield stress of highly concentrated MR fluids.

It is well-known that the yield stress of MR fluids scales as 0.5 1.5 0 s MH  , with the magnetic field intensity H, for applied fields of medium intensity (far from the initial linear region of the magnetization curve of the material and from the saturation region) [START_REF] Bossis | Ferrofluids[END_REF]. Here 0  is the vacuum permeability and s M the saturation magnetization of the magnetic particles of the MR fluid. We used this scaling and, as observed, the curves for the different applied field strengths approximately superimpose in a single master curve (Figure 8b). For this purpose, we used the saturation magnetization of the Fe-CC particles, s M = 870.9 ± 1.4 kA m -1 [14]. We can thus conclude that in concentrated MR fluids, for each imposed value of the shear rate, the resulting shear stress approximately scales with H 1.5 .

IV. CONCLUSIONS

We have shown that in the pre-yield regime a steady-state flow is impossible in highly concentrated suspensions of magnetic particles subjected to an applied magnetic field. At a given value of the imposed shear rate within the pre-yield regime, instead of a quasi-constant value of the shear stress at large enough times, we obtained shear stress oscillations of well-defined amplitude and period. Rheograms obtained from the experimental measurements displayed an N-like shape, with the initial (at very low shear rate) increasing region of stress and the decreasing one, corresponding to the values of shear rates for which oscillations were observed. The final increasing region of stress, at high enough values of the shear rate, corresponded to the regime of developed flow, for which no oscillations were obtained. Both stress oscillations and the N-like shape disappeared in the absence of applied field and were not observed under any experimental condition for suspensions of moderately low particle concentration.

Concerning the effect of the measuring system gap thickness on the obtained Nlike rheograms, there was overlap of the rheograms in the developed flow regime.

However, in the pre-yield regime, the shear stress tended to increase as the gap thickness was decreased. This increment was very likely related to an enhancement of particle structures under compression, as supported by the increase of normal force that we obtained experimentally. Finally, we have also shown that for a given value of the imposed shear rate the resulting shear stress increased as 1.5 H , with the magnetic field intensity H, for the range of magnetic fields used in our work.

Figure 1 .

 1 Figure 1. Sketch of the rough (serrated) geometry used with MCR 300 (Physica-Anton

Figure 2 .

 2 Figure 2. Time evolution of the shear stress for different imposed values of the shear

Figure 3 .

 3 Figure 3. Shear rate and shear stress plotted as a function of shearing time for MR

Figure 4 .

 4 Figure 4. Shear rate and shear stress plotted as a function of shearing time for a MR

Figure 5 .

 5 Figure 5. Rheograms of MR suspensions. (a) In the absence of applied magnetic field,

Figure 6 .

 6 Figure 6. Time evolution of the normal force (a) and the shear stress (b), as the gap

Figure 7 .

 7 Figure 7. Time evolution of the normal force as the gap thickness between the rotating

Figure 8 .

 8 Figure 8. Rheograms of a suspension containing 50 vol.% of iron particles under the

ACKNOWLEDGEMENTS

This work was supported by project FIS2013-41821-R, MINECO, Spain; by project 3.12.2014/К, Program of Ministry of Science and Education of the Russian Federation; by the Act 211 Government of the Russian Federation, contract № 02.A03.21.0006; by grants of RFFI 13-02-91052, 13-01-96047 (Ural) and 14-08-00283; and by project PICS 6102 CNRS/Ural Federal University. L.R.-A acknowledges financial support from University of Granada.

Under an applied magnetic field, the rheograms of suspensions containing 50 vol.% of iron particles presented a marked N-like shape, characterized by an initial increase of the shear stress with the shear rate, followed by a decrease of the shear stress as the shear rate was further increased, and a final, rather strong, increase of the shear stress for values of the shear rate higher than approx. 10 s -1 (Figure 5b). Similar N-like rheograms were also obtained for measuring plates with rough surfaces, as well as for measuring plates with a magnetic layer coating (Figure 5c). However, as observed, values of the shear stress for a given value of the shear rate were considerably higher for measuring plates with a magnetic layer coating. This enhancement of the MR effect can be explained as a result of the better transmission of forces to the working plate, in agreement with previous works [START_REF] Laun | Wall material and roughness effects on transmittable shear stresses of magnetorheological fluids in plate-plate magnetorheometry[END_REF][START_REF] Jonkkari | Effect of the plate surface characteristics and gap height on yield stresses of a magnetorheological fluid[END_REF]. On the other hand, the N-like shape was not so significant in the absence of applied magnetic field (Figure 5a). Furthermore, it was not appreciably detected in rheograms of suspensions containing 10 vol.% of iron under an applied magnetic field (Figure 5d).

Note that similar N-like rheograms to these of suspensions of 50 vol.% of iron particles under applied magnetic field were reported for laponite suspensions [START_REF] Pignon | Thixotropic colloidal suspensions and flow curves with minimum: identification of flow regimes and rheometric consequences[END_REF]. Very recently, we also observed similar rheograms for suspensions of magnetic nanofibers [START_REF] Zubarev | N-Like rheograms of suspensions of magnetic nanofibers[END_REF], although in this previous work the decreasing region of shear stress was obtained at much higher values of the shear rate (around 10 s -1 ). The decrease of the shear stress with the shear rate in the case of magnetic nanofibers was connected to the existence of strong friction between particles.

Effect of the gap thickness of the measuring system


Laurent Daudet
email: daudet@lam.jussieu.fr

Stéphane Molla
email: molla@cmi.univ-mrs.fr

Bruno Torrésani

B Torresani
email: torresan@cmi.univ-mrs.fr

Towards a Hybrid Audio Coder

come

Introduction

Digital audiophonic signal coding has become an important issue in many application areas. Among the most popular approaches, transform coding has received particular attention during the recent decades, as the rapid development of hardware and the discovery of novel mathematical approximation techniques has made it particularly efficient. Transform coding starts with an expansion into a suitably chosen orthonormal basis of the spaces of signals. In the (functional) approximation stage, only a (small) subset of the coefficients is retained, and encoded (after quantization, which is the second point where approximation comes into play).

Although much may be done at the level of quantization, we focus here on functional approximation, which may be thought in two different ways:

• Linear approximation: in this approach, the retained coefficients correspond to a fixed subset of a suitably chosen basis a . • Non-linear approximation: the retained basis functions and coefficients are selected adaptively, so that the corresponding truncated expansion minimizes a given distorsion (some norm of the approximation error). Non-linear approximation schemes automatically outperforms the linear schemes (using the same basis and the same number of retained coefficients) in terms of distorsion, but introduces an extra costs in terms of encoding: the retained set of basis functions being not fixed a priori, the corresponding information (the addresses of retained coefficients) has to be encoded as well.

To cope with the problem of address encoding for non-linear approximation, the concept of structured approximation has to be introduced: it stems from the fact that for given classes of signal (and/or functional spaces), and accordingly chosen orthonormal bases, the significant coefficients have a natural tendancy to cluster around some given types of structures (lines/tubes, trees,...) in their index space. Exploiting such information yields substantial gains in the addresses coding, and thus in the performances of encoders. In addition, this often improves the quality of signal modeling, as we shall see on more specific examples below.

The coding schemes we shall outline in the present paper are based upon structured non-linear approximations, with an additional non-linearity, introduced for the following reason. Audio signals (like other signal classes, images [START_REF] Meyer | Multilayered Image Representation: Application to Image Compression[END_REF] ,...) may be thought as compound objects, containing significantly different features [START_REF] Serra | A System for Sound Analysis/Transformation/Synthesis based on a Deterministic plus Stochastic Decomposition[END_REF] : mainly tonals (usually termed partials in the audio and speech literature), transients, and additional residual components. The adequate orthonormal bases for transform coding these components (at least the first two) are significantly different, and we therefore model the signal as a linear sum of two different components, each being transform coded with an adequate basis. To our knowledge, our approach is fundamentally new, in that it does not rely on a prior segmentation of the signal into different components [START_REF] Levine | Audio Representations for Data Compression and Compressed Domain Processing[END_REF][START_REF] Verma | A Perceptually Based Audio Signal Model With Application to Scalable Audio Compression[END_REF] : superimposition is prefered to segmentation. Preliminary a One generally uses an orthonormal basis which approximates well the Karhunen-Loève basis of the signal (modeled as a random signal), under some additional constraints (essentially, the existence of efficient algorithms) results validating this approach have already been presented in previous publications [START_REF] Daudet | Hybrid models for audio signals encoding[END_REF] . We shall discuss here different levels of hybrid modeling, outline important problems (theoretical and practical) encountered in this context, and present estimation schemes developed so far.

Approximations and structures

The significantly different components present in audiophonic signals may generally be parsimoniously represented using suitable bases. This remark has led several authors to propose to expand such signals on dictionaries, obtained as unions of these bases. These dictionaries generally form quite redundant systems, and raise the problem of finding the optimal signal expansion on such a system (among all such expansions), optimality being understood in terms of parsimony. Focusing on the particular application to audio signals, and limiting ourselves to transient and tonal features, we are naturally led to consider a generic redundant dictionary made out of two orthonormal bases [START_REF] Vetterli | Wavelets and subband coding[END_REF][START_REF] Wickerhauser | Adapted Wavelet Analysis from Theory to Software[END_REF] (typically a wavelet basis {ψ n , n = 0, . . . N -1} and an MDCT basis {w m , m = 0, . . . N -1}, and signal expansions of the form

x = λ∈Λ α λ ψ λ + δ∈∆ β δ w δ + r , (1)
where Λ and ∆ are (small, this being the sparsity assumption) subsets of the index sets, hereafter termed significance maps. In addition, we also introduce a residual signal r, which is not sparse with respect to the two considered bases (a spread residual), and is to be neglected or described differently.

Assuming that such an expansion has been obtained, the corresponding tonal layer and transient layer of the signal x are defined by

x ton = δ∈∆ β δ w δ , x tran = λ∈Λ α λ ψ λ . (2)
Several approaches have been proposed to perform such estimations, see for example [START_REF] Donoho | Uncertainty principles and ideal atomic decompositions[END_REF][START_REF] Elad | A generalized uncertainty principle and sparse representations[END_REF][START_REF] Gribonval | Sparse representations in union of bases[END_REF] but most of these are not necessarily adapted when it comes to practical implementation in a coding perspective. On one hand, it is not clear that the corresponding algorithms are compatible with practical constraints, in terms of CPU and memory requirements. Also, models exploiting solely sparsity arguments cannot capture one of the main features of these signal classes, namely the persistence property: significant coefficients have a tendency to form "structured sets": lines or tubes of MDCT coefficients for tonal signals, and trees 17 of wavelet coefficients for transient signals.

We outline below three approaches of increasing complexity, that illustrate the benefits of introducing such "structural information" into the approximation techniques.

Exploiting solely sparsity: N-term approximation

We first describe briefly a first approach exploiting only sparsity arguments [START_REF] Berger | Removing noise from music using local trigonometric bases and wavelet packets[END_REF][START_REF] Daudet | Hybrid models for audio signals encoding[END_REF] . The tonal layer x ton in ( 2) is essentially obtained by thresholding of the MDCT transform (after suitable weighting,that corrects for the overall decay of MDCT coefficients 5 ). The transient layer x tran is obtained similarly, by thresholding the wavelet transform of the nontonal signal x nton = x -x ton . The values of the thresholds, and therefore the corresponding significance maps ∆ and Λ, are determined by the required bit rate. In order to allocate appropriate amount of bits to the two layers, the relative proportion may be pre-estimated for each time frame, using appropriate techniques [START_REF] Molla | Determining local transientness of audio signals[END_REF] . The resulting algorithm reads as follows:

Algorithm 1: Within each time frame:

(1) Pre-estimate the relative importance of the tonal and transient layers, and corresponding bit rates.

(2) Expand the signal on an MDCT basis, and pick the largest coefficients (as estimated in Step 1). Reconstruct the tonal signal x ton and the nontonal part x nton = x -x ton .

(3) Expand x nton on a wavelet basis, and select the largest coefficients (according to the rule given in Step 1). This generates the transient layer x tran , and the residual signal x res = x nton -x tran .

Structured N-term approximation

As stressed in the introduction of this paper, the performances of the coding scheme are greatly improved if "structure" information is implemented into the encoder [START_REF] Daudet | Hybrid models for audio signals encoding[END_REF] . During the estimation of the tonal layer, significant MDCT coefficients are retained only if they satisfy some time-persistence property. Let β δ = β k,ν be the MDCT coefficients, where k relates to time (index of the transform window w k ), and ν relates to frequency. Sine waves around frequency bin ν are fully characterized by the coefficients β kν , β k,ν-1 and β k,ν+1 , and the local Fourier spectrum can be well approximated [START_REF] Daudet | MDCT analysis of sinusoids: exact results and applications to coding artifacts reduction[END_REF] , with a smoothed MDCT spectrum on a given window w t :

βk,ν = β 2 k,ν + (β k,ν+1 -β k,ν-1 ) 2 1 2 (3)
Now, a width-3 tube T k1,k2 ν = {β k,ν-1 , β k,ν , β k,ν+1 } k=k1...k2 is retained if its averaged (over time) smoothed spectrum :

σ p [T k1,k2 ν ] = 1 |T k1,k2 ν | k=k1..k2 w(ν)| βk,ν | p (4)
exceeds some fixed minimum value. Here, |T k1,k2 ν | = k 2 -k 1 + 1 denotes the time duration of the tube, w(ν) are frequency-dependent weights (e.g. related with the absolute threshold of hearing), and p is a constant selecting different types of tubes (strong short vs. long weak).

After estimation and substraction of the tonal layer, the transient layer is estimated in a similar way: the nontonal signal is expanded into a wavelet basis, and only significant wavelet coefficients satisfying some scale persistence property are retained. As is well known, wavelet coefficients are naturally organized into dyadic trees. A branch B of the wavelet coefficient tree is retained if the following modulus of regularity

κ q,s [B] = 1 |B| λ∈B 2 j(λ)s |α λ | q , (5)
exceeds a given maximum value (here, |B| denotes the length of the branch B, and j(λ) denotes the scale parameter corresponding to the index λ).

The constants s, q characterize the considered type of transients, as they weight coefficients corresponding to different scales.

Algorithm 2: Within each time frame (1) Pre-estimate the relative importance of the tonal and transient layers, and corresponding bit rates.

(2) Expand the signal on an MDCT basis, and estimate the tubes. Reconstruct the tonal layer x ton and the non-tonal part x nton = x -x ton .

(3) Expand x nton on a wavelet basis, and estimate the branches of the tree. Reconstruct the transient layer x tran , and the residual x res = x nton -x tran .

Hybrid structured hidden Markov model

The algorithms outlined above do not rely on any modeling for the considered signals, and it is therefore difficult to make any a priori estimate regarding their performances. For that purpose, a model of hybrid waveform audio model has been developed 15 , which implements hidden Markov chains (resp. trees) of MDCT (resp. wavelet) coefficients. To be more specific, a tonal signal is still modeled as in (2), and the (significant) coefficients β δ , δ ∈ ∆ are modeled as N (0, σ2 δ ) independent random variables. The index δ = (k, ν) is a pair of time-frequency indices and the significance map ∆ is characterized by a "fixed frequency" Markov chain involving two hidden states, hence by a set of initial frequencies ν 1 , . . . ν N and transitions matrices P1 , . . . PN (one for each frequency bin). Therefore, we model the MDCT coefficients as a mixture model of two gaussians processes, characterizing heavy-tailed like processes. Globally, the tonal model is characterized by the set of matrices Pn , and the variances σ 2 δ of the two states, which are assumed to be time invariant, and on which additional constraints may be imposed.

A similar model, using Hidden Markov trees of wavelet coefficients 3 may be developed to describe the transient layer in the signal, as in (2). To model the scale persistence of large wavelet coefficients of transients, the significant wavelet coefficients {α λ , λ ∈ Λ} of the signal are modeled as N (0, σ 2 λ ) random variables. The index λ = (j, k) is a pair of scaletime indices and the significance map Λ is characterized by a "fixed time frame" Markov chain, hence by "scale to scale" transition matrices P j (with additional constraints ensuring that significant coefficients inherit a tree structure, see below.) In addition, we impose that the transition from a "non-transient" state towards a transient state is impossible, so that non significant coefficients will not be followed by significant ones. Then, a transient structure is defined within each frame as a rooted binary tree of significant coefficients.

The transient model is therefore characterized by the variances of wavelet coefficients in Λ and Λ c , and the persistence probability, for which estimators may be constructed. The transient states estimation itself is also performed via classical methods (involving well known Expectation-Maximization technique).

Numerical examples

We now illustrate the several approaches described above on numerical examples, which we use to emphasize the differences between them. For the sake of simplicity, we limit ourselves to a single sound example, namely a small piece extracted from a song of Ben Harper. The signal is exhibited at the top of the plots in Figures 4, 5 and6. In all three situations, the overall number of retained coefficients (tonal + transient) was set to 5% of the number of samples of the original signal, therefore a functional compression ratio of 20. As may be seen, in all cases, the tonal layer may be seen as a copy of the original signal with a somewhat "smoothed" envelope, while the transient layer has essentially captured the "attacks". However, there are some noticeable differences between the three results, on which we elaborate below. The main difference comes from the difference between N-term and structured N-term approximations. We have represented in Figures 4 and5 nal (top), tonal layer (middle) and transient layer (bottom). The original signal features both tonals ("horizontal structures") and transients ("vertical features"), and it may be seen than these are significantly different in Figures 4 and5. As expected, the structured N-term approximation yields more "persistent" structures than the simple N-term approximation. This leads to a better resolution of the tonal and transient features. Concerning the signal and layers themselves, we remark that the signal at hand is mostly tonal, as confirmed by the transientness index plotted in Figure 3. Therefore, most of the energy of the signal is contained in the tonal layer, and the three methods do not make much difference regarding that particular component. Overall, between 80% and 85% of coefficients were spend on the tonal layer. While the tonal layers appear quite similar at first sight, the transient layers do exhibit significant differences. We notice that the simple N-term approximation yields transient signal present at all time, while the structured methods (structured N-term and Markov) yield more lacunary transients, which seem more relevant for the signal at hand.

The structure of the transient layer obtained in the structured N-term method appear fairly terser than in the other situations. This is presumably due to the estimation procedure, which is in that case quite simple, and does not really exploit the tree structure: let us recall that in order to keep a simple algorithm, we have limited the detection of transient structures to branches rather than complete trees. In this respect, an exploration of trees should probably be prefered at this stage, yielding higher computational costs. This point is currently under study. At the opposite, the Markov A closer examimation of the residual signal also shows that the results obtained using the Markov model are more satisfactory: indeed, the residual looks "more stationary" in that situation than in the other two cases, which is good in a coding perspective [START_REF] Daudet | Hybrid models for audio signals encoding[END_REF][START_REF] Molla | An Hybrid Audio Scheme using Hidden Markov Models of Waveforms[END_REF] . Again, this is probably partly due to a poor estimation of the transient layer in the N-term and structured N-term approaches.

Conclusions

We have briefly outlined in this paper three approaches of increasing complexity for hybrid encoding of audio signals. The results presented here show that hybrid modeling and structured approximation are definitely suitable techniques for audio coding. We have outlined the main differences between the three approaches, and their practical consequences. Among other conclusions, it appears clearly that a good modeling of wavelet significance trees is necessary for good estimation of the transient layer. Let us also point out that the pre-estimation of tonal and transient bit rates is also an important ingredient, and that our current approach probably needs further refinement. More complete results in a signal encoding context will be presented in a forthcoming publication. Additional illustrations and sounds available at http://www.cmi.univ-mrs.fr/~torresan/papers/Chongqing

Algorithm 3 :

 3 Within each time frame (1) Pre-estimate the relative importance of the tonal and transient layers, and the corresponding bit rates. (2) Learning stage: expand the signal on an MDCT basis, and estimate the best parameters for the model. Estimation stage: estimate relevant lines. Reconstruct the tonal layer x ton and the non-tonal part x nton = x -x ton . (3) Learning stage: Expand x nton on a wavelet basis, and estimate the corresponding parameters of the Markov model. Estimation stage: estimate transient trees. Reconstruct the transient layer x tran , and the residual x res = x nton -x tran .

Figure 1 .

 1 Figure 1. MDCT modulus of original signal (top), tonal (middle) and transient layers for the Ben Harper audio signal based on N-term approximation.

Figure 2 .

 2 Figure 2. MDCT modulus of original signal (top), tonal (middle) and transient layers for the Ben Harper audio signal based on structured N-term approximation.

Figure 3 .

 3 Figure 3. Transientness index for the Ben Harper signal

Figure 4 .

 4 Figure 4. Hybrid expansion of the Ben Harper audio signal based on N-term approximation; from top to bottom: original, tonal layer, transient layer, residual

Figure 5 .

 5 Figure 5. Hybrid expansion of the Ben Harper audio signal based on structured N-term approximation; from top to bottom: original, tonal layer, transient layer, residual

14

 14

Figure 6 .

 6 Figure 6. Hybrid expansion of the Ben Harper audio signal based on Hybrid hidden Markov model; from top to bottom: original, tonal layer, transient layer, residual

* This work is supported by the French Ministry of Research (contract "ACI Jeunes Chercheurs"). † Work partially supported by hassip.


Jason Chevrie

Alexandre Krupa

Marie Babel

Needle Steering Fusing Direct Base Manipulation and Tip-based Control

Works on robotic needle steering often consider either tip-based control of flexible beveled tip needles or control of the bending of symmetric tip needles. In this paper a control law for needle steering which uses both direct manipulation of the needle base and the control of the tip deflection, is proposed. A 3D model of a beveled tip needle is used to determine the base motion that is needed to obtain the desired tip motion. The needle-tissue interaction model uses local virtual springs placed along the needle shaft and takes into account the interaction forces at the bevel. Online estimation and update of this model parameters is performed via visual feedback. The low level controller uses the task function framework to allow control of the tip velocity. Additionally duty cycling method is used if a reduction of the natural deflection of the needle tip is needed. Finally an experimental targeting task in a gelatin phantom is presented for preliminary validation of the framework. The method provided sub-millimeter accuracy on a target that would be unreachable using only tip-based control.

I. INTRODUCTION

Surgical interventions using needles have been increasingly used over the last decades. One well-known advantage of needles is the reduced trauma that they cause in the tissue. This allows non-invasive operations, thus reducing the recovery time needed after the operation. However, accurate needle steering in soft tissue remains a challenging task for the clinicians. Indeed, this surgical act is sensitive to any deformations of the tissues, target displacements or needle bending. At the same time the clinician has reduced visibility since no direct vision is possible. Naturally, the effect of mistargeting can be traumatic. It is typically the case in brachytherapy, where misplacement of the radioactive seeds can cause the destruction of healthy tissues. Robot assisted procedures can thus be a great help in this task and have been an active research field for many years [START_REF] Reed | Robot-assisted needle steering[END_REF].

In this context, the bending of the needle was at first considered a drawback that had to be reduced. Indeed, the steering of a flexible needle requires from the clinician an even more accurate comprehension of the interaction between the needle and the tissues [START_REF] Reed | Robot-assisted needle steering[END_REF]. However, it has then been shown that using flexible needles could allow target reaching that would otherwise be unreachable by a straight line [START_REF] Webster | Nonholonomic modeling of needle steering[END_REF]. A robotic assistance of this challenging task is thus even more necessary.

Pioneer work was conducted by DiMaio et al. [START_REF] Dimaio | Needle steering and motion planning in soft tissues[END_REF] on flexible needle steering in soft tissues. A finite element model was 1 Jason Chevrie is with Université de Rennes 1 and IRISA, France, jason.chevrie@irisa.fr 2 Alexandre Krupa is with Inria Rennes -Bretagne Atlantique and IRISA, France, alexandre.krupa@inria.fr 3 Marie Babel is with Insa Rennes and IRISA, France, marie.babel@irisa.fr used to model both needle and tissues. The needle motion and tissue deformations could be predicted when the needle was moved by its base. Nevertheless, real time modeling could not be achieved because of the model complexity. A simpler framework was proposed by Glozman et al. [START_REF] Glozman | Image-guided robotic flexible needle steering[END_REF] [START_REF] Neubach | Ultrasound-guided robot for flexible needle steering[END_REF]. The needle-tissue interaction was there modeled only locally with virtual springs placed along the needle shaft. In both methods the needle tip trajectory is controlled by using the transverse motion of the base to deform the tissues and bend the needle. This provides a good steerability at the beginning of the insertion (low depth in tissue). However it can result in large tissue stress and possible tissue damage once the needle is deeply inserted.

From the observation that a beveled needle tends to bend when inserted [START_REF] Okamura | Force modeling for needle insertion into soft tissue[END_REF], most recent works mainly focus on the steering of very flexible needles using this property. It has been shown by Webster et al. [START_REF] Webster | Nonholonomic modeling of needle steering[END_REF] that the trajectory of the needle tip could be modeled by simple kinematic models, such as the unicycle or bicycle ones. Many control schemes have then been proposed using these models. In [START_REF] Rucker | Sliding mode control of steerable needles[END_REF] sliding mode control is proposed such that the bevel is always oriented toward the target. Hauser et al. [START_REF] Hauser | Feedback control for steering needles through 3d deformable tissue using helical paths[END_REF] use the helical trajectories obtained when simultaneously inserting and rotating the needle. When the rotation speed is large with respect to the insertion speed, straight trajectories can be obtained. The duty cycling control method [START_REF] Minhas | Modeling of needle steering via duty-cycled spinning[END_REF] allows a linear control of the effective curvature of the trajectories. This facilitates the reduction of the complexity of some high level planners, such as the Rapidly exploring Random Tree planner [START_REF] Patil | Needle steering in 3-d via rapid replanning[END_REF]. The method was used with an online estimation of the natural curvature for a targeting task in [START_REF] Moreira | Needle steering in biological tissue using ultrasound-based online curvature estimation[END_REF].

The basic assumption when using the kinematics models is that the needle shaft exactly follows the tip trajectory. This is the case when using a very flexible needle embedded in a hard medium. In practice, however, the tissue can be rather soft and the needles commonly used by clinicians not completely flexible. Misra et al. [START_REF] Misra | Mechanics of flexible needles robotically steered through soft tissue[END_REF] proposed a complete mechanics based model to predict the tip deflection when the needle is inserted. Contrary to the kinematics models, this model takes into account the tissue deformations due to the bending of the needle inside the tissue. In [START_REF] Abayazid | Integrating deflection models and image feedback for real-time flexible needle steering[END_REF] this model was simplified by adding the same virtual springs approach as in [START_REF] Glozman | Image-guided robotic flexible needle steering[END_REF] and used to control the needle trajectory.

Our first contribution in this paper is the design of a model that simulates the evolution of a beveled tip needle in 3D space and in real time. The needle tissue interaction is modeled by using an extension of the virtual springs model [START_REF] Glozman | Image-guided robotic flexible needle steering[END_REF] to 3D. A similar bevel model to the one proposed in [START_REF] Misra | Mechanics of flexible needles robotically steered through soft tissue[END_REF] is added to take into account the force generated at the needle tip during insertion. The model we proposed is then used to derive the relationship between the motion of the tip and the motion of the base of the needle.

The second contribution is a new control law based on the proposed needle model. It allows the steering of the needle tip toward a desired direction. The originality of our approach is to use transverse motion control to increase the reachable region and to switch to a duty cycling method when a reduction of the natural curvature is needed, i.e when rather straight trajectories must be achieved. For the rotation along the shaft, we use sliding mode control [START_REF] Rucker | Sliding mode control of steerable needles[END_REF] to orient the bevel in such a way that the needle tends to bend in the desired direction. We also added a secondary task to reduce the tissue deformations when possible.

The paper is organized as follows: our 3D model of the needle is presented in section II. In section III we present the different parts of our control law. Our experimental setup and results are then shown in section IV. Finally, conclusions and future works are shown and discussed in section V.

II. NEEDLE MODELING WITH VIRTUAL SPRINGS

In this section we introduce the model that we use to estimate the behavior of the needle when the motion of the base is controlled. The model proposed by Glozman et al. [START_REF] Glozman | Image-guided robotic flexible needle steering[END_REF] is first extended: modifications are made in such a way that high bending angles of a 3D needle can be modeled. This allows the modeling of highly flexible needles trajectories in space. The tip model proposed by Misra et al. [START_REF] Misra | Mechanics of flexible needles robotically steered through soft tissue[END_REF] is also combined to consider the case of a beveled tip.

A. Spring Model

The interaction between the needle and the tissue is modeled locally with N virtual springs placed all along the needle shaft. Each spring i is defined by a rest position p 0,i , a stiffness K i and a plane P i that contains p 0,i (see Fig. 1a). The rest position of the spring corresponds to the rest location of the tissue if no needle is pushing on it. Each spring exerts a force on the needle at the location where the needle crosses the plane P i . We use here the model in a quasi-static case. So we assume that the frictional part, that is tangential to the needle, can be neglected. The springs are only used to model the normal forces F s,i applied on the needle shaft. For this, the orientation of the planes P i are constantly updated such that the exerted force stays normal to the needle. The force can thus be expressed according to

F s,i = K i (p n,i -p 0,i ), (1)
with p n,i the point of the needle which crosses the plane P i .

We consider here that the tissue stiffness per unit length K T around the needle is constant all along the needle. The stiffness of each spring can thus be adapted in such a way that the equivalent tissue stiffness per unit length remains the same all along the needle, that is

K i = K T l i (2)
where l i is the length of the needle that is supported by the spring K i .

B. Bevel tip

As proposed by Misra et al. in [START_REF] Misra | Mechanics of flexible needles robotically steered through soft tissue[END_REF], we choose here to model the efforts exerted at the tip with triangular loads distributed on each side of the tip. Let α be the bevel angle, β the cut angle, b the bevel length, a the length of the bottom edge of the needle and O the center of the needle located just before the bevel (see Fig. 1b). The normal force F tip and moment M tip exerted at the point O can be expressed as

F tip = K T a 2 2 tan β - K T b 2 2 tan(α -β) cos α, (3)
M tip = - K T a 3 6 tan β + K T b 2 2 tan(α -β) a 3 cos(α) - b 6 sin(α) 2 .
Knowing the orientation of the base frame F b (see Fig. 1a), the orientation of the tip frame F t around the tip axis can be deduced. We choose here to assume that the torsional bending of the needle can be neglected. This way, we can obtain the orientation of F t by integrating the local rotation around the out of plane vector all along the needle.

C. Needle Model

The needle is modeled by a succession of N + 1 segments such that the extremities of the segments lay on the planes of the virtual springs, except for the needle base that is fixed and the needle tip that is free. Each segment is approximated in 3D using a polynomial curve S i (l) of order n so that

S i (l) = M i [ 1 l . . . l n ] T , (4)
where i ∈ [1, N + 1] is the segment index and S i (l) ∈ R 3 is the position of a point of the segment at the curvilinear coordinate l ∈ [0, L i ] with L i the total length of the segment.

M i ∈ R 3×(n+1
) is a matrix containing the coefficients of the polynomial curve. In our case, real-time performance is expected: so the model complexity should be low. A third order polynomial is thus chosen, since it is the lowest sufficient order for which the mechanical equations can directly be solved.

To maintain adequate continuity properties of the needle, second order continuity constraints are added, namely defined as

S i (L i ) = S i+1 (0), (5)
dS i dl l=Li = dS i+1 dl l=0 , ( 6
)
d 2 S i dl 2 l=Li = d 2 S i+1 dl 2 l=0 . ( 7
)
The total normal force F i at the extremity of the segment i can be calculated from the sum of the forces exerted by the springs located from this extremity to the needle tip, so that where Π i stands for the projection onto the plane P i . The projection is used to remove the tangential part of the force and to keep only the normal part. This normal force introduces a constant shear force all along the segment, such that EI

F i = Π i   F tip y tip + N j=i F s,j   , (8)
d 3 S i dl 3 (l) = -F i , (9)
with E the needle Young's modulus and I its second moment of area. For a hollow circular needle, I can be calculated from the outer and inner diameter [START_REF] Glozman | Image-guided robotic flexible needle steering[END_REF], d out and d in respectively, with the well-known formula

I = π 64 (d 4 out -d 4 in ). (10)
Finally the moment due to the bevel force gives the following boundary condition:

EI d 2 S N +1 dl 2 l=L N +1 = -M tip y tip . (11)
The above modeling provides a total of 12 × (N + 1) unknown variables and 12 × (N + 1) equations. Like in [START_REF] Glozman | Image-guided robotic flexible needle steering[END_REF], it is thus possible to solve this problem as a simple linear problem.

III. HYBRID CONTROL FOR FLEXIBLE NEEDLE

In this section, we propose a control law that couples both base manipulation and tip based controls for the steering of the needle. At the best of our knowledge it is the first time that a method fusing these controls is proposed. The control enables the steering of the needle tip in a desired direction, while trying to reduce the tissue deformation. When the needle is aligned with the target, duty cycling allows to straighten the trajectory and reduce the needle deflection.

Let v d t be the desired translation velocity (dim(v d t ) = 3) of the needle tip expressed in the tip frame. In practice this velocity will be given by another algorithm that generates the desired trajectory.

A. Direct Base Manipulation

Let V b = v b ω b
T be the velocity screw vector of the base, with v b = v b,x v b,y v b,z T being the velocity vector and ω b = ω b,x ω b,y ω b,z T the angular velocity vector. Similarly we define the velocity screw vector of the tip, V t = v t ω t T . Using the previously defined model, we numerically derive a manipulation Jacobian t J b . It represents a local approximation of the relation between the base and tip screw vectors given by

V t = t J b V b . ( 12
)
The translational velocity of the tip is thus given by

v t = t J b,v V b , (13)
where t J b,v is the upper half part of t J b .

In order to reduce the tissue deformation while obtaining the desired tip trajectory, we propose to use the redundancy formalism of the task function framework [START_REF] Siciliano | A general framework for managing multiple tasks in highly redundant robotic systems[END_REF]. This allows considering many different tasks that have to be fulfilled by the system. These tasks are hierarchically ordered and the control law is computed in such a way that low priority tasks do not induce any perturbations on higher priority tasks.

In our case, we define three tasks. A first task e 1 of dimension 3 is used to control the trajectory of the needle tip. This task has the highest priority and is defined such that ė1 = v t . The desired behavior for this task is here ẽ1 = v d t . A second task e 2 of dimension 1 is used for the orientation of the bevel and has a medium priority. If we do not orient the bevel, additional needle bending can be necessary to compensate for the natural tip deflection. This would increase the tissue stress and is of course not desirable. To naturally reduce the tissue deformation and needle bending that is required to steer the tip in the desired direction, the bevel should face the opposite direction. In this way, the tip axis naturally orients itself toward the desired velocity vector when the needle is simply inserted. Thereby, the second task is defined as e 2 = σ, where σ is the angle between the current bevel orientation and the orientation given by v d t (see Fig. 1c), i.e.

σ = - π 2 -atan2(v d t,y , v d t,x ). ( 14
)
We define the desired behavior of the task as ẽ2 = ω b , where ω b is given by a smoothed sliding mode controller [START_REF] Rucker | Sliding mode control of steerable needles[END_REF]. This control law allows to steer σ toward zero and is such as:

ω b = -sign(σ) ω z,max if |σ| ≥ π 4 -4 π σ ω z,max if |σ| < π 4 ( 15
)
where ω z,max is the maximum rotation speed allowed. Finally a third task e 3 of dimension 1 is used for active reduction of the tissue deformations. We give this task the lowest priority and define it as e 3 = d m , where d m is the mean tissue displacement. This is estimated using the model proposed in section II by averaging the elongation lengths of the virtual springs

d m = 1 L ins N i=1 l i p n,i -p 0,i , (16)
where L ins is the total length of the needle that is inserted in the tissue. To obtain a smooth reduction of the tissue deformation, we choose the desired behavior of the task as ẽ3 = -λd m , where λ is a coefficient that influences the convergence speed. Since the optimal directions of these three tasks may be in opposition (mainly the first and third tasks), we use a singularity-robust formulation for the control law [15][16]. This way, the final control law for the needle base motion is given by

V b = t J + b,v v d t + P 1 t J + b,ωz ω b -λP 12 J + dm d m , (17)
where the + symbol stands for the pseudo-inverse operator, t J b,ωz is the sixth line of the manipulation Jacobian and J dm is the Jacobian of the mean tissue deformation with respect to the base screw vector. P 1 is here the orthogonal projection matrix operator onto the null space of t J b,v and P 12 is the projection matrix operator onto the null space of the augmented Jacobian J 12 , namely:

J 12 = t J b,v t J b,ωz , (18)
P 1 = I 6 -t J + b,v t J b,v , (19)
P 12 = I 6 -J + 12 J 12 , (20)
where I 6 is the 6 × 6 identity matrix.

B. Duty Cycling Control

The duty cycling method is used to reduce the natural curvature obtained when inserting the needle. This has been extended to 3D using different approaches in [START_REF] Wood | Algorithm for threedimensional control of needle steering via duty-cycled rotation[END_REF] and [START_REF] Krupa | A new duty-cycling approach for 3d needle steering allowing the use of the classical visual servoing framework for targeting tasks[END_REF]. The ratio between the natural curvature K nat of the needle and the effective curvature K ef f has been shown in [START_REF] Minhas | Modeling of needle steering via duty-cycled spinning[END_REF] to be linearly dependent on the duty cycle DC, i.e. the ratio between the two phases durations:

K ef f = K nat (1 -DC) (21)
On a kinematics point of view, DC actually corresponds to the ratio between the insertion length of the phases. In our case, we define a fixed insertion length L DC for the cycles and allow the use of different insertion speeds for the two phases. This way we can avoid the infinite cycle duration obtained when DC goes to zero as stated in [START_REF] Krupa | A new duty-cycling approach for 3d needle steering allowing the use of the classical visual servoing framework for targeting tasks[END_REF]. We use duty cycle control when the desired tip velocity v d t is almost aligned with the needle tip axis. This is detected by comparing the angle θ d between v d t and the tip axis (Fig. 1c) with the angle θ DC that would be obtained during a cycle with duty cycle equal to zero, i.e. with an insertion phase only. θ d and θ DC are expressed as

θ d = atan2 v d t,x 2 + v d t,y 2 , v d t,z , (22)
θ DC = K nat L DC . ( 23
)
If θ d < θ DC , it means that the needle would overshoot the current desired direction in less than a cycle length. In that case it is better to reduce the effective curvature of the needle such that it aligns with the desired direction, i.e with

K ef f = θ d L DC . ( 24
)
Using the same approach as in [START_REF] Krupa | A new duty-cycling approach for 3d needle steering allowing the use of the classical visual servoing framework for targeting tasks[END_REF], we also define an accumulation angle for the rotation around the tip axis equal to σ (see eq. ( 14)), such that the bevel is in the adequate direction before starting the translation phase.

Finally, since the needle can be bent, the optimal insertion direction of the base V b,ins may be different from the needle base axis. Here we choose to compute it as the base motion that maximizes the tip insertion in the axial direction, V b,ins = t J + b,vz , where t J b,vz is the third line of the manipulation Jacobian.

IV. EXPERIMENTAL RESULTS

In this section we present the experimental setup used to validate the control method and the results.

A. Experimental setup

A needle is fixed to the end effector of a 6 DOFs Adept Viper S650 robot. We use a stainless steel biopsy needle with internal and external diameters of 0.48 mm and 0.7 mm respectively. The length of the needle that can bend, i.e. the part that is outside the needle holder, is 12.6 cm long. The needle Young's modulus is assumed to be that of steel, i.e. 200 GPa. The insertion is done in a homemade gelatin phantom embedded in a transparent plastic container. For a validation purpose, the visual feedback is here obtained using two optical cameras. The needle is detected and tracked in real time by an image processing algorithm that we briefly introduce in the following. A white screen monitor is used to provide some light in the background and enables the observation of the needle through the gelatin. Figure 2e shows a picture of this experimental setup. The intrinsic parameters of both cameras were calibrated as well as the extrinsic parameters between the cameras. In this way, stereo-vision is used to retrieve the 3D position of the needle observed in both camera images. The accuracy of this tracking system, corresponding to the size of the pixels, is around 0.25 mm. Before the experiments, the gelatin container was removed and the shape of the needle was detected. As the position of the needle is known with respect to the robot, the transformation between the robot frame and the stereo-vision system was estimated. Note that we plan in the future to replace this visual system by a 3D ultrasound probe, thus providing directly a 3D position without the need for stereo-vision calibration.

B. Model update from visual feedback

The pose of the needle base in the model is updated using the odometry of the robot and the fact that the orientation of the bevel with respect to the robot is known at the beginning of the procedure. The needle is tracked in each image and the 3D shape of the needle that is visible in both images is computed. As the needle is clearly visible in the images, we use here a simple tracking algorithm that consists in estimating the needle shape with a third order polynomial defined by four control points. The position of the control points is updated and filtered with a Kalman filter at each image acquisition such that they follow the needle shape and its tip. We don't describe the algorithm further here because it is not the subject of this paper and it should be adapted to a more relevant imaging modality in the future.

To have a good modeling of the tissue, virtual springs must be regularly added all along the modeled needle shaft. We choose here to add a new spring when the length of the last needle segment reaches a fixed threshold of 2.5 mm. When adding a spring at the needle tip, we define the rest position of the spring at the location of the tip. This way the rest position represents the initial location of the tissue when it is just being cut by the bevel. To fit the model to the reality, this position is obtained from the real needle tip using the visual feedback.

The model also requires some parameters as inputs. The parameters related to the needle, such as its geometry and Young's modulus, are already known and should not vary significantly during the insertion procedure. The tissue stiffness per unit length K T , however, must be estimated. We propose here to do an online estimation of this parameter. The best value is computed at each image acquisition using a least square approach such that the shape of the model needle fits the measured needle shape.

C. Experiments

Four experiments were performed to validate our method. The stiffness per unit length of the model was initialized in each case with K T =10000 N/m 2 . At the beginning of each experiment, the needle is placed such that it is normal to the surface of the gelatin and its tip slightly touches it. The insertion point is shifted between the experiments in a way that the needle can not cross a previous insertion path. The needle is first inserted along a distance of few millimeters in the gelatin to allow the manual initialization of the tracking algorithm in the images. Then the insertion procedure is started and is stopped when the target is not in front of the needle tip anymore.

In a first experiment, a virtual target is defined before the beginning of the insertion such that it is aligned with the needle and placed at a distance of 8 cm from the tip. An insertion is then performed without activating the controller at an insertion speed of 1 mm/s, i.e. the needle is simply inserted along its shaft direction. The insertion is stopped when the needle tip passes the target location. Fig. 2a shows the view of the front camera at the end of the experiment. We can see that the target is missed due to the natural deflection of the needle. Fig. 3a shows the 3D lateral distance between the needle tip axis and the target. It can be seen that the needle passes almost 8 mm away from the target at the end of the insertion.

In a second experiment, the designed controller is activated. The input velocity vector v d t for the controller is chosen such that it is always pointing toward the target and has a norm of 1 mm/s:

v d t = 0.001 p tar p tar , (25)
with p tar the position of the target in the tip frame, deduced from the visual feedback. As can be seen on Fig. 3c the controller allows to reach the target with a sub-millimeter accuracy, thus validating the framework. Two other experiments were performed similarly to the previous ones, except that the target is shifted 1 cm away from the initial tip axis. For the insertion without the controller (third experiment), the target is placed in the bevel direction (opposite of the y axis on Fig. 1b) such that the needle naturally bends toward the target when inserted. For the last experiment, the controller is activated and the target is placed to forms an angle of 135 degrees with the bevel direction to be sure that the controller has to rotate the needle to align the bevel with the target.

We can see on figure Fig. 2b that the target can not be reached when the controller is not active. Fig. 3b shows that the 3D lateral distance between the target and the tip axis tends to decrease but the target is finally missed by 5 mm. A pure insertion actually corresponds to a saturated duty cycling control with duty cycle DC equal to 0. Using only the natural tip deflection to control the needle is thus not enough for this targeting task because of the low curvature of 3. Measure of the 3D lateral distance between the needle tip axis and the target. Note that the measure is noisier at the beginning, since errors on the measure of the tip axis orientation are amplified when the target is far from the tip. the needle. Since our controller uses additional base motions, it allows to reach the target as can be seen on Fig. 2c and2d. Fig. 3d shows that the target could also be reached with sub-millimeter accuracy. This way the designed controller allows to increase the reachable space compared to the sole use of tip-based control methods. Table I shows a summary of the final lateral targeting error between the tip and the target. Note that the axial error (in the needle direction) was under 0.2 mm for each experiment, which corresponds to the accuracy of the vision system.

V. CONCLUSIONS AND FUTURE WORK

In this paper we proposed a new control method for flexible needle steering that combines direct base manipulation and tip based control methods. The direct base manipulation control is generated thanks to the use of a 3D model of a beveled tip needle that gives the adequate motion of the needle base to obtain a given motion of the tip. An update of the model is performed online using visual feedback. The duty cycling method is jointly used to reduce the needle natural curvature and to obtain relatively straight trajectories when needed. We observed that our controller allows to obtain a good targeting accuracy and is able to reach targets that are not reachable using only tip-based control.

As a next step, the cameras will be replaced by 3D ultrasound probe to fit with the clinical context. The tracking algorithm should be adapted for ultrasound images, as for example by using the one presented in [START_REF] Chatelain | 3d ultrasound-guided robotic steering of a flexible needle via visual servoing[END_REF]. One drawback with the current formulation of the control law is that little deformation reduction is achieved when the task comes in opposition with the targeting task. Least squares optimization could be used to find a better control, but at the expense of greater computation load. We also plan to extend the update of the model to take into account the displacements of the tissues that are due to the physiological motion of the patient.

  Fig. 1. Needle modeling

Fig. 2 .

 2 Fig. 2. Final view of the front camera at the end of the experiments. The crosses represent the target. Refer to the accompanying video for real time view of the insertions. (e) Experimental setup

  Fig.3. Measure of the 3D lateral distance between the needle tip axis and the target. Note that the measure is noisier at the beginning, since errors on the measure of the tip axis orientation are amplified when the target is far from the tip.

TABLE I FINAL

 I LATERAL ERROR in MM

		No control Active control
	Aligned target	7.6	0.6
	Shifted target	4.9	0.1