M Heinisch

T Voglhuber-Brunnmaier

E K Reichel

I Dufour

B Jakoby

Electromagnetically Driven Torsional Resonators for Viscosity and Mass Density Sensing Applications

Keywords: Torsional, resonator, cylinder, shear wave, viscosity, mass density

In this contribution a conceptual study for torsional oscillators, which are electromagnetically driven and read out, is presented. The aim is to experimentally investigate the basic feasibility of a torsional resonator with application to viscosity and mass density sensing in liquids. Such a device is particularly interesting as cylindrical, torsional resonators for fluid sensing applications are hardly reported but unlike many other devices, yield pure shear wave excitation in the liquid. The design of first conceptual demonstrators for measurements in air as well as in liquids and their benefits and disadvantages are discussed in detail. A closed form as well as a reduced order model and measurement results obtained with first demonstrators are presented.

Introduction

Recently, we investigated various resonant sensors for liquid viscosity and mass density, see e.g., [START_REF] Jakoby | Miniaturized sensors for the viscosity and density of liquids -performance and issues[END_REF], which were particularly designed to be operated in the low kilohertz range. Amongst these devices, in-plane oscillating platelets, emitting mainly shear waves into the sample liquids, were investigated e.g., in [START_REF] Reichel | Analysis and experimental verification of a metallic suspended plate resonator for viscosity sensing[END_REF][START_REF] Abdallah | Measurement error estimation and quality factor improvement of an electrodynamic-acoustic resonator sensor for viscosity measurement[END_REF] where millimeter sized metal platelets were used. In [START_REF] Riesch | A suspended plate viscosity sensor featuring inplane vibration and piezoresistive readout[END_REF][START_REF] Cerimovic | Sensing viscosity and density of glycerol-water mixtures utilizing a suspended plate mems resonator[END_REF] similar miniaturized devices have been implemented in silicon technology. Generally used shear oscillating resonators, such as shear oscillating quartz crystals [START_REF] Martin | Characterization of a quartz crystal microbalance with simultaneous mass and liquid loading[END_REF] and the aforementioned in-plane oscillating platelets, have in common that, besides the desired shear waves, compressional waves are also radiated into the liquid. These pressure waves result e.g., from non-uniform shear displacement [START_REF] Beigelbeck | A two-dimensional analysis of spurious compressional wave excitation by thickness-shear-mode resonators[END_REF][START_REF] Voglhuber-Brunnmaier | Efficient spectral domain formulation of loading effects in acoustic sensors[END_REF], the resonators' finite thicknesses, and spurious out-of-plane modes e.g. from the plate itself or of supporting beams. Potential candidates for resonators which only emit shear waves into the test fluid are cylindrical torsional oscillators. First, such pure shear wave emitting devices are of special interest from a rheological point of view, when it comes to the analysis of complex liquids such as viscoelastic liquids. For pure shear wave deformation, the liquid can be described by a complex-valued shear modulus or a complex valued viscosity, which, in general, are frequency dependent quantities [START_REF] Macosko | Rheology, Principles, Measurements and Applications[END_REF]. With oscillatory measurements, the values obtained at a number of discrete frequencies can be used to obtain a rheological spectrum. A second merit for cylindrical, torsional resonators is that, due to their geometry, they resemble concentric cylinder rheometers. Thus, they yield comparable measurements but extend the measurable frequency range of conventional rheometers which is usually limited at approximately 100 Hz. Third, cylindrical, torsional oscillators allow experimental comparison of measurement results obtained with the above mentioned in-plane oscillating resonators, which excite spurious compressional waves. Thus, by means of such a comparison the impact of these compressional waves could be experimentally estimated.

In-plane oscillating devices, oscillating rotational disks [START_REF] Rahafrooz | Characterization of rotational mode disk resonator quality factors in liquid, Frequency Control and the European Frequency and Time Forum (FCS)[END_REF] and other resonators for viscosity and mass density sensors such as cantilevers [START_REF] Dufour | The microcantilever: A versatile tool for measuring the rheological properties of complex fluids[END_REF][START_REF] Vančura | Analysis of resonating microcantilevers operating in a viscous liquid environment[END_REF], quartz tuning forks [START_REF] Toledo | Application of quartz tuning forks and extensional microresonators for viscosity and density measurements in oil/fuel mixtures[END_REF] and vibrating bridges [START_REF] Riesch | A novel sensor system for liquid properties based on a micromachined beam and a low-cost optical readout[END_REF][START_REF] Etchart | Mems sensors for density-viscosity sensing in a low-flow microfluidic environment[END_REF] are based on a similar operational principle. Usually, the devices' frequency responses, containing a characteristic resonant mode, are recorded upon immersion in a sample liquid. The change of evaluated resonance frequencies and quality factors are then related to the liquid's mass density and viscosity.

In this work a feasibility study for electromagnetically driven and read out torsional oscillators for operation in liquids is presented. Two concepts for actuation and readout discussed in Sec. 2 were manufactured and the experimental analysis in air allowed investigating the benefits and disadvantages of both approaches. The analysis of measurements obtained with three different torsional spring diameters and various spring lengths allowed designing a demonstrator for measurements in liquids, which is explained in detail in Sec. [START_REF] Abdallah | Measurement error estimation and quality factor improvement of an electrodynamic-acoustic resonator sensor for viscosity measurement[END_REF]. A closed form model for the resonator, relating the input to the output voltage, is derived in Sec. 4. In Sec. 5 evaluated results of measurements obtained with the cylindrical torsional resonator in ten different liquids are shown. Furthermore, the sensor's sensitivity and cross sensitivity to temperature are evaluated and compared to other sensor concepts.

Torsional resonator

Concept

The basic idea of the torsional resonator is to excite a cylinder to torsional vibrations by means of Lorentz forces acting on sinusoidal currents in a constant external magnetic field. For recording the device's frequency response, the excitation current's frequency is swept over a frequency range containing the resonant fundamental mode and simultaneously measuring the motion induced voltage on an electrical conductor, following the torsional oscillation. To implement this idea, a rotor (bobbin) is mounted on torsional springs, where two different principles were realized and compared. In the first approach, the bobbin carries one single coil which is used for both, excitation and read out, where in the second approach, two separate coils are used for these tasks.

Conceptual investigation in air

To investigate the functional principle of the electromagnetic torsional oscillator and to estimate the achievable range of resonance frequencies and signal strengths, a single coil type bobbin has been mounted and investigated on three torsional springs at various spring lengths. For this, a setup as depicted in Figs. 1(a) and 1(b) has been used. There, one hundred turns of a 80 µm thick copper wire were wound on a 3D printed bobbin with 8 mm in diameter and 22 mm in length, which was attached to tungsten rods with diameters of 0.58 mm, 1.6 mm and 2 mm serving as torsional springs. For each torsional spring, the same rotor was used so that only the effect of different torsional spring lengths and diameters could be examined. For attaching, the rotor was affixed with screws to the torsional springs, which were rigidly clamped at their ends with fibreglass blocks through which the torsional spring lengths could be adjusted. The bobbin was placed in a magnetic field (denoted with B) provided by neodymium permanent magnets and set to torsional oscillations by means of Lorentz forces on sinusoidal currents (I in ) in the coil which was connected to a signal generator (V g and 50 Ω output resistance) and a series resistance R s = 100 Ω which was used to limit the excitation current, to prevent from non-linear deflections. By sweeping the excitation current's frequency, the oscillator's frequency response can be recorded. In this case, for the sake of straight forward manufacturing, the coil's ends of the 80 µm thick copper wire were kept long enough for direct connection with the excitation and readout electronics. However, for a stable resonator, this wiring approach is not adequate. In the demonstrator used for measurements in liquids which will be presented in Sec. 3 this drawback has been overcome by connecting the coil's ends to the torsional springs for electrical connection of the resonator.

Figure 2 shows the results of recorded frequency responses and evaluated resonance frequencies in air for tungsten rods with diameters of 0.58 mm and 1.6 mm for various spring lengths in comparison with theoretical results also presented in [START_REF] Heinisch | Concept study on an electrodynamically driven and read-out torsional oscillator[END_REF]. These results proved the basic feasibility of the concept and allowed designing a torsional resonator in the desired frequency range.

Single versus double coil setup

The benefits and disadvantages of using single or double coil setups were experimentally investigated. The cross sec- tions of the rotors as well as the electrical equivalent circuits for both cases are depicted in Figs. 1(b) and 1(c). For the first case, a single coil with 100 turns, a length of 20 mm and a width of 8 mm is used for both, excitation and read out and for the second case, two coils with 50 turns each are used separately for these purposes. The voltage of the signal generator V g was 0.1 V and the series resistance R s was 100 Ω.

For the single coil setup, the rotor is excited to torsional vibrations by means of Lorentz forces on sinusoidal currents I in . The oscillation of the coil in the presence of the external magnetic field causes a motion induced voltage V M on the excitation coil itself. This motion induced voltage is proportional to the the velocity of the oscillation and thus an appropriate quantity for measuring the oscillation. However, the measurable output voltage V out consists of an additional voltage resulting from the coil's impedance related voltage drop. The coil's impedance is modeled as a serial connection of a resistance R c and an inductance L c . For the manufactured rotor R c = 21.3 Ω and L c = 215 µH were obtained by fitting the parameters in the recorded frequency response from 100 Hz to 100 kHz. The output voltage was measured with a lock-in amplifier (with an input resistance of 10 MΩ) and reads

V out (ω) = V M (ω) + (R c + j ω L c )I in (1)
in complex notation where j = √ -1, ω is the angular frequency and the time dependence e j ω t is suppressed.

To reduce this impedance-related offset voltage, a second coil is used in the double coil setup for readout. This second coil follows the motion of the oscillating cylinder and thus, as in the single coil setup, a voltage is induced. However, in this case due to electrical crosstalk resulting from the current I in in the excitation coil, an additional voltage V ct is induced in the measuring coil. Assuming that the output current I out is negligible, the output voltage reads

V out (ω) = V M (ω) + j ω M I in ( 2
)
where M is the mutual inductance describing the inductive coupling from the excitation to the readout coil. For the rotor used in this experiment M = 480 µH was obtained. Both types of rotors have been manufactured and mounted on a tungsten rod with 0.58 mm diameter and 5 cm length for each torsional spring. The recorded frequency responses covering a frequency range of 100 Hz to 10 kHz as well as detailed plots of the fundamental resonance are shown in Fig. 3 for both cases. The resonance peak is more than twice as high in case of the single coil setup compared with the double coil setup. (12.4 mV in contrast to 5.5 mV for the double coil setup.) This results from the fact that for the single coil setup, the number of coil turns is twice as high which yields a higher effective driving force in the rotor and a proportional to the coil turns higher induced voltage. Thus, assuming the same mechanical conditions for both setups, the voltage peak of the single coil setup • is that the output voltage is subjected to a slight crosstalk voltage at these frequencies. This crosstalk voltage increases proportionally to the excitation frequency as it can be observed in Fig. 3(b) which is also considered in Eq. 2.

Demonstrator for measurements in liquids

The use of a second coil did not show a major advantage in the obtained signals compared to the setup where only one coil is used for excitation and readout. However, using only one coil has a significant advantage for electrical connection. In this case, the electrical conductive torsional springs can be used for this task. For the double coil setup which requires the connection of four wires, the usage of tubes would be necessary (serving as torsional springs) through which the coils' wires have to be threaded which might be a rather tedious manufacturing process. Fig. 4 shows drawings and photographs of the torsional resonator with which first measurements in ten different liquids were obtained.

As for the single coil type rotor, which was investigated in air, one hundred turns of a copper wire with a diameter of 80 µm were wound on a 3D printed bobbin. Two tungsten rods with a diameter of 0.58 mm serve as torsional springs with a spring length of 14 mm each. Each rod is brazed to a brass cylinder with 10 mm diameter and 10 mm height. Both tungsten rods carry a brazed brass platelet and are put into the bobbin on both sides ensuring no contact of both torsional springs. Each end of the coil is glued with electrical conductive glue to one of these brass platelets which are furthermore used for a form-fitted connection with epoxy resin poured into both sides of the bobbin which in turn is put into a POM-C (polyacetal) plastics jacket with an outer diameter of 10 mm and 25 mm length. As mentioned before, it is essential that the tungsten rods do not get in contact in order not to shortcut the coil. The torsional resonator is put into a milled POM-C plastics frame with a well for containing the sample liquids. The brass cylinders are clamped with brass screws, which in turn are used for further electrical connection for power supply and read-out. Two neodymium magnets are used for providing a magnetic field necessary for excitation and read-out based on Lorentz-forces.

Modeling

Mechanical model 4.1.1. Structural mechanics

Considering the rotor's moment of inertia J 0 , the torque resulting from the tungsten rods' torsional stiffnesses M T , the excitation torque M ex , a torque M c , representing the intrinsic dissipative losses of the resonator and a torque M F accounting for the torque due to the liquid loading, the principle of momentum equilibrium yields, cf. Fig. 5 J 0 d 2 ϕ(t)

dt 2 + 2 M T + M c + M F = M ex ( 3
)
where ϕ is the twisting angle of the rotor. For a cylindrical rod with shear modulus G, radius r s , length l s and twisting angle ϕ the torsional torque reads [START_REF] Weaver | Vibration Problems in Engineering[END_REF]

M T = k * T ϕ. ( 4
)
l s l r M T J M T M c M F M ex
is the torsional spring constant of one torsional spring and k T = 2 k * T is the spring constant of both torsional springs in our case. Dissipative losses are considered by a loss coefficient c 0 and the associated loss torque is

M c = c 0 dϕ dt . (6)

Fluid dynamics

As we are facing a rotational problem, the solution of the fluid forces acting on the oscillating cylinder contains Hankel functions which make the obtained equations difficult to interpret. For this reason, the fluid forces acting on the cylinder are approximated by one dimensional shear waves of an in-plane oscillating plate which yields a negligible error if the cylinder's radius r c is significantly larger than the penetration depth [START_REF] Landau | Fluid Mechanics[END_REF], also termed decay length [START_REF] Martin | Characterization of a quartz crystal microbalance with simultaneous mass and liquid loading[END_REF]

δ = 2η ρ ω (7)
where η and ρ are the liquid's viscosity and mass density, respectively and ω is the angular frequency of the oscillation.

A comparison of the relative deviations ε of the solutions for the cylindrical and the planar case over the ratio r c /δ are depicted in Fig. 6. The complete derivation for the fluid forces is only given for the planar case in this contribution. For aqueous liquids at 1 kHz the penetration depth is 18 µm approximately. Considering the cylinder's radius of the manufactured demonstrator r c = 5 mm it follows r c /δ = 280 and thus ε = 3.8 • 10 -3 which substantiates the applicability and the validity of the approximated solution.

Under the assumptions of an infinitely extended plate oscillating only in x-direction (thus imposing only shear stress in x-direction), that gravitational forces are negligible, and that the liquid is incompressible, the equation of motion and the shear stress T xz in liquids can be expressed as a one dimensional problem and reads [START_REF] Martin | Characterization of a quartz crystal microbalance with simultaneous mass and liquid loading[END_REF][START_REF] Heinisch | A suspended plate in-plane resonator for rheological measurements at tunable frequencies[END_REF]:

ρ ∂ 2 u x ∂t 2 = ∂T xz ∂z with T xz = η ∂ 2 u x ∂z ∂t . (8)
Here, ρ is the mass density of the liquid, u x is the displacement of the liquid in x-direction, t is the time variable. Substituting T xz in the equation of motion, transforming the problem to the frequency domain assuming a time dependence e j ω t and solving the linear differential equation of second order yields the solution for the x-displacement propagating in z-direction

u x (z, t) = u 0 e - z δ e -j z δ -ω t (9)
where u 0 is the amplitude of the oscillation. With this solution for u x , the shear stress at the liquid-solid interface (i.e., z = 0) in the frequency domain reads on an oscillating cylinder and an in-plane oscillating plate, respectively. For both cases, the same surface interacting with the liquid is considered. In case of the solution for the cylinder, Hankel functions are obtained which make an intuitive interpretation of the obtained equations difficult. For high ratios of the cylinder's radius r c and the penetration depth δ, the difference between both solutions is negligibly small.

T xz (z = 0) = u x (1 -j) η ρ ω 3 2 . ( 10
)
Expressing the planar deflection in rotational form, i.e. u x → ϕ r c and assuming M F = -r c T A c , where A c is the cylinder's surface, the torque acting on the oscillating cylinder can be expressed as

M F = -ω 2 η ρ 2 ω r 2 c A c + jω η ρ ω 2 r 2 c A c ϕ(ω). (11)
Eqs. 4 and 6 are substituted in Eq. 3. and after transformation of the obtained expression to the frequency domain, Eq. 11 can be considered in the equation of the moment equilibrium

-ω 2 J + k T + jω c ϕ(ω) = M ex ( 12
)
where J and c have been introduced and read

J = J 0 + η ρ 2 ω r 2 c A c and c = c 0 + η ρ ω 2 r 2 c A c .( 13
)
As the torsional resonator is operated close to its (fundamental) angular resonance frequency ω 0 , these expressions can be approximated by

J ≈ J 0 + η ρ 2 ω 0 r 2 c A c and c ≈ c 0 + η ρ ω 0 2 r 2 c A c .( 14
)
Substituting these relations into Eq. 12 yields a closed form model of the oscillating cylinder in liquids and is beneficial for design purposes. However, as the model consists of a relatively 4.2. Reduced order models 4.2.1. One dimensional shear waves Combining factors and variables in Eq. 14 in single factors and neglecting the ω 0 dependence allows writing J and c in reduced form as follows

J ≈ J 0 + J * ηρ √ η ρ and c ≈ c 0 + c * ηρ √ η ρ. ( 15
)
With this, the angular speed Ω(ω) = jωϕ(ω) can be expressed as

Ω(ω) = M ex /c 1 + j Q ω ω 0 -ω 0 ω ( 16
)
with angular resonance frequency

ω 0 = k T J = 1 J 0k + J * ηρk √ η ρ (17)
and quality factor

Q = √ J k T c = 1 ω 0 k T c = J 0k + J * ηρk √ η ρ c 0k + c * ηρk √ η ρ . ( 18
)
where the subscript k denotes that the factors in Eq. 15 have been divided by k T to obtain these solutions for ω 0 and Q.

For resonators, which are readout via a motion induced voltage (which is proportional to the oscillator's velocity, see Sec. 4.3.2) the maximum peak frequency and the frequency of free, undamped oscillations (usually termed ω 0 ) are identical. For this reason, we may call ω 0 resonance frequency. For read-out principles evaluating the resonator's deflection using e.g. optical methods [START_REF] Baù | Cost-effective system for the characterization of microstructures vibrating in out-of-plane modes[END_REF], the frequency of free, undamped oscillations and the maximum peak frequency of deflection amplitude are not identical.

Generalized model

In [START_REF] Heinisch | Reduced order models for resonant viscosity and mass density sensors[END_REF] a similar reduced order model was presented which considers not only one dimensional shear waves but convexly shaped, oscillating objects in general and reads

ω 0 = 1 m 0k + m ρk ρ + m ηρk η ρ ω 0 (19)
and

Q = 1 ω 0 • 1 c 0k + c ηk η + c ηρk √ ω 0 η ρ , (20)
where m xk and c xk are coefficients and m ρk as well as c ηk are zero in case of pure one-dimensional shear waves.

Simplified generalized equations

Equation 19 is an implicit equation for ω 0 , which makes an exact evaluation of ω 0 for given η and ρ difficult. Therefore, the equations are simplified by considering that the frequency dependence of certain parameters occurring in the analysis is negligible as they are virtually constant within the bandwidth of the resonant system. Doing so, we obtain the following expressions for the angular resonance frequency

ω 0 = 1 m 0k + m ρk ρ + m * ηρk √ η ρ (21)
and the quality factor

Q = m 0k + m ρk ρ + m * ηρk √ η ρ c 0k + c ηk η + c * ηρk √ η ρ . ( 22
)
Thus if ω 0 and Q for given η and ρ using Eqs. 19 and 20 have to be calculated, numerical, (e.g., iterative) methods can be used. Alternatively, the above simplified expressions give remedy. Conversely, if η and ρ have to be determined from measured f r = ω 0 /(2 π) and Q, Eqs. 19 and 20 can be used directly.

The advantage of such reduced order models is that in comparison with completely closed form models, only a few variables have to be determined by means of a parameter fit.

In Sec. 5 these models are fit to experimental results. By means of a comparison of the fitted models it will be shown that the generalized models fit the measured data better than the model considering shear waves only. This supposedly is due to the fact that in addition to one dimensional shear waves, spurious effects such as liquid trapping [START_REF] Martin | Sensing liquid properties with shear-mode resonator sensors[END_REF][START_REF] Jakoby | Physical sensors for liquid properties[END_REF] may be present.

Electrical model 4.3.1. Excitation torque

The Lorentz force on a charge q moving with velocity v in the presence of a magnetic field with magnetic flux density B and an electric field E is (bold variables denote vectors) [START_REF] Griffiths | Introduction to Electrodynamics[END_REF]

F L = q (v × B + E). ( 23
)
Assuming perpendicular conditions between the excitation current I ex , i.e. the charges moving in a wire of length l r , and the external magnetic field, it follows that the force acting on the wire is given by

F ex = B l r I ex . ( 24
)
With this, the excitation torque can be expressed as

M ex ≈ 2 N B r r l r I ex ( 25
)
where N is the number of turns of the excitation coil and r r is the radius of the coil, see Fig. 7.

Motion induced voltage

Using Faraday's law of induction and the equation for Lorentz forces, the motion induced voltage on the coil can be expressed as:

V M (t) = - coil (v × B) • ds. ( 26
)
Assuming the magnetic flux density to be symmetric in the air gap, constant along the coil's length and aligned in parallel with the coil for ϕ = 0, see Fig. 7, it follows

V M (t) = 2 N r r cos ϕ dϕ dt B l r . (27)
For small twisting angles, the motion induced voltage in the frequency domain using the angular speed Eq. 16 reads

V M (ω) ≈ 2 N r r B l r Ω(ω). ( 28
)
This equation for the motion induced voltage is used for modeling the measurable output voltage V out for the single as well as for the double coil setup in Eqs. 1 and 2.

Experimental investigation

Measurements in liquids

The demonstrator discussed in Sec. 3 has been used to experimentally test the principle operation of such a torsional resonator in liquids and to investigate the effect of different mass densities and viscosities on recorded frequency responses. Amplitudes as well as evaluated quality factors and resonance frequencies for two different sets of liquid series are depicted in Fig. 8. There, the results of the fitted models performing a linear fitting procedure, see [START_REF] Chong | An introduction to optimization[END_REF], for f r = ω 0 /(2 π) and Q using Eqs. 17, 18, 21 and, 22 are depicted as well. The values for the examined liquids' viscosities η and mass densities ρ (which were determined with an Anton Paar SVM 3000) as well as the mean values of the associated, evaluated resonance frequencies f r and quality factors Q are given in Tab. 1. The SVM 3000 features a reproducibility of 0.35 % for viscosity and 0.0005 g/cm 3 for mass density. The measurements were performed at 25 • C. The values for f r and Q were evaluated from the recorded frequency responses using a fitting algorithm described in [START_REF] Niedermayer | Methods for the robust measurement of the resonant frequency and quality factor of significantly damped resonating devices[END_REF]. The used liquid series are the socalled viscosity series and mass density series in the following. The first is a mixture of isopropanol and acetone covering a viscosity range of 0.21 mPa•s to 2.05 mPa•s for mass densities of roughly 0.78 g/cm 3 . The percentage of mass isopropanol m I and mass acetone m A is given in the left column of Tab. 1. The mass density series cover a range of mass densities of 0.79 g/cm 3 to 1.01 g/cm 3 for viscosities of roughly 1 mPa•s and were prepared using acetone, isopropanol, ethanol, DI-water and glycerol. After mixing, the liquids were investigated with the SVM 3000. These two liquid series are used to separately investigate experimentally the effect of varying viscosities or mass densities on f r and Q.

The results obtained with both liquid series clearly show the effect of mass density and viscosity on the experimentally obtained values for f r and Q. Deviations of the measured from the modeled values as well as slight variations in f r and Q can be observed. These supposedly originate from imperfections in the prototype setup. These imperfections are associated with the fact that first, during cleaning and refilling of the sensor, the A c c e p t e d M a n u s c r i p t resonator may have been detuned due to an insufficiently stable clamping, second, the open experimental well, which allows evaporation of the liquid during the measurements, and third, geometrical imperfections such as e.g., surface roughnesses and a not perfect alignment of the oscillating cylinder with the torsional springs. For Q, both models yield approximately the same results, where for f r , the generalized model clearly yields better results. This finding can be explained by surface roughness, which yields so-called liquid trapping [START_REF] Martin | Sensing liquid properties with shear-mode resonator sensors[END_REF][START_REF] Jakoby | Physical sensors for liquid properties[END_REF] and thus a higher sensitivity to mass density. This effect might be beneficial, as the effect of η and ρ on f r and Q cannot be separated with devices yielding pure one-dimensional shear waves without liquid trapping. Despite the discussed drawbacks in the presented setup, a rough estimation of the device's sensitivity can be made. We define the relative sensitivity of a quantity X(y i ) to an independent variable y i as

Viscosity Series T = 25 • C m I m A +m I / % η/(mPa•s) ρ/(g/cm 3 ) f r /Hz Q η/(mPa•s) ∆ η/(mPa•s) ∆ η rel ρ/(g/cm 3 ) ∆ ρ/(g/cm
S X,y i = ∂X ∂y i • y i X (29)
where X either stands for f r or Q and y i for η and ρ, yielding four different sensitivities which are evaluated from the fitted generalized model Eqs. 21 and 22. The four sensitivities S f r ,η , S f r ,ρ , S Q,η and S Q,ρ are not constant but depend on η as well as on ρ and thus are depicted as bands in Fig. 9 for the experimentally evaluated range of viscosities and mass densities. The upper boundary of the bands in the plots on the left hand side (i.e., S f r ,η and S Q,η ) are the values for ρ = 1.01 g/cm 3 , the lower boundary for ρ = 0.79 g/cm 3 . On the right hand side (i.e., S f r ,ρ and S Q,ρ ) the upper and lower boundaries of the bands represent the evaluated values for η = 0.21 mPa•s and η = 2.05 mPa•s, respectively. In Fig. 9 the sensitivities obtained with the torsional resonator are depicted in comparison with the sensitivities achieved with circular and rectangular cross sectioned steel tuning forks oscillating at 400 Hz in liquids approximately, a U-shaped wire sensor ( f r ≈ 930 Hz) [START_REF] Heinisch | A u-shaped wire for viscosity and mass density sensing[END_REF], and a quartz tuning fork ( f r ≈ 32.7 kHz) [START_REF] Toledo | Application of quartz tuning forks and extensional microresonators for viscosity and density measurements in oil/fuel mixtures[END_REF].

Measurement accuracy

Equations 19 and 20 were used to calculate η and ρ for evaluated f r and Q. The fitted model parameters using a linear fitting procedure described in [START_REF] Heinisch | Reduced order models for resonant viscosity and mass density sensors[END_REF][START_REF] Chong | An introduction to optimization[END_REF] are given in Tab. more accurate measurement results for viscosity. However, the purpose of this work was a feasibility study of torsional, resonant viscosity and mass density sensors, their modeling and first designs of promising demonstrators. In [START_REF] Heinisch | Application of resonant steel tuning forks with circular and rectangular cross sections for precise mass density and viscosity measurements[END_REF] the advantages and capability of resonant viscosity and mass density sensors using conventional steel tuning forks with circular and rectangular cross sections was demonstrated.

Cross sensitivity to temperature

Mechanical resonators can show a significant dependence of their resonance frequency to temperature. This dependence is mainly due to the thermal expansion of the resonator and the temperature dependence of the resonator's Young's modulus. For doubly clamped structures such e.g. bridges [START_REF] Etchart | Mems sensors for density-viscosity sensing in a low-flow microfluidic environment[END_REF][START_REF] Riesch | Micromachined Viscosity Sensors[END_REF] and straight wire resonators [START_REF] Seibt | Schwingdrahtviskosimeter mit integriertem ein-senkkörperdichtemessverfahren für untersuchungen an gasen in größeren temperatur-und druckbereichen[END_REF][START_REF] Heinisch | Tunable resonators in the low khz range for viscosity sensing[END_REF], the cross sensitivity of the resonance frequency to temperature becomes large if significant thermal pre-stresses are induced. The dependence of the resonance frequency to temperature can be positive as well as negative. In [START_REF] Larsen | Ultrasensitive string-based temperature sensors[END_REF], micro-machined vibrating strings are used as ultrasensitive temperature sensors. For the case of resonant viscosity sensors, this cross sensitivity directly limits the sensor's accuracy and thus should be kept as low as possible. The dependence of the torsional resonator's resonance frequency has been evaluated in a temperature range from 5 • C to 60 • C in 5 • C temperature steps. Figure 10 shows the magnitude of measured frequency responses as well as evaluated resonance frequencies at these temperatures. In Tab. 3 a comparison of the resonance frequency's dependence to temperature is given for different resonators. The first five sensors, i.e a U-shaped tungsten wire sensor, U (W), [START_REF] Heinisch | A u-shaped wire for viscosity and mass density sensing[END_REF], a U-shaped gold coated silicon res- onator, U(Si, Au), [START_REF] Stifter | Mems heterodyne amf detection with capacitive sensing[END_REF], a steel tuning fork, TF (steel), a silicon cantilever, CL (Si), [START_REF] Wasisto | Airborne engineered nanoparticle mass sensor based on a silicon resonant cantilever[END_REF] and a gold coated silicon cantilever, CL (Si, Au), [START_REF] Sandberg | Temperature and pressure dependence of resonance in multi-layer microcantilevers[END_REF] are singly clamped structures and thus show a relatively small dependence of the resonance frequency to temperature. The last two sensors i.e., the torsional resonator, TR, and a straight tungsten wire sensor, SW (W), [START_REF] Heinisch | Tunable resonators in the low khz range for viscosity sensing[END_REF] are doubly clamped structures and thus, show a much higher dependence of the resonance frequency to temperature. However, this dependence is much lower for the torsional resonator compared to the straight wire sensor. To reduce this cross sensitivity, using only one torsional spring instead of two may be a possible approach.

Sensor Table 3: Cross sensitivities of the resonance frequency (fundamental mode) to temperature. f 0 : nominal resonance frequency, ∆ f r = f r (T 0 )f r (T 1 ): resonance frequency shift, ∆T = T 0 -T 1 difference of temperatures T 0 and T 1 , ∆ f r,rel = ∆ f r / f 0 : relative resonance frequency shift, 't.w.' designates this work.

f 0 ∆ f r /∆T ∆ f

Conclusion and outlook

A first study of the electromagnetically driven and read out torsional oscillator in air showed the advantages of using one single coil for both, excitation and readout. The advantages are mainly based on the manufactural benefit that in case of using one coil (instead of two) only two wires have to be electrically connected for power supply and readout which in this case is implemented by contacting the torsional springs. The feasibility of such a torsional resonator for viscosity and mass density sensing applications has been shown by measurements in ten different liquids. For design purposes, a complete analytical model comprising the structural and fluid mechanics and relating output to input signals has been elaborated. Therefrom, reduced order models have been derived which relate resonance frequency and quality factor to viscosity and mass density, respectively. The advantage of such models is that only four to six factors have to be determined to describe the sensors characteristics. The comparison of both fitted model results shows that the generalized model yields significantly better results for f r which can be explained by liquid trapping due to surface roughnesses of the cylinder shell. A first estimation of the device's sensitivity and its cross sensitivity of its resonance frequency has been made and compared to other resonant viscosity and mass density sensors.

Regarding future work, a more reliable clamping and housing has to be designed and manufactured aiming to obtain more stable and accurate results. By means of further experiments, the maximum measurable viscosities and resolution have to be investigated.

Page 10 of 12

Figure 1

 1 Figure 1: a) Principle and photographs of first demonstrators allowing to record frequency responses in air, b) cross section and electrical equivalent circuit of the demonstrator with a rotor carrying a single coil for excitation and readout by means of Lorentz-forces. I in : input (excitation) current, B: external magnetic field, V g : voltage of the signal generator, 50 Ω: output resistance of the signal generator, R s : series resistance, R c , L c : coil's resistance and inductance, V M : motion induced voltage, 10 MΩ: input resistance c) cross section and electrical equivalent circuit of the rotor carrying two coils. One coil is used for excitation, the other for read out., V ct : induced voltage due to inductive crosstalk.

Figure 2 :

 2 Figure2: Examples for recorded frequency responses in air using a single coil for excitation and read-out. These experiments were performed to estimate achievable signal strengths and resonance frequencies for different spring diameters and spring lengths.

Figure 3 :

 3 Figure 3: Frequency responses measured with the single coil (a) and the double coil (b) setup. The recorded spectra between 100 Hz and 10 kHz as well as detailed plots of the resonance are shown for both setups.

Figure 4 :

 4 Figure 4: a) Exploded view of the torsional resonator used for measurements in liquids. The coil is electrically connected via the torsional springs and the brazed brass cylinders, which are clamped with brass screws to which connecting wires are soldered. b) shows an illustration of the assembled, clamped setup, the direction of the external magnetic field as well as the torsional oscillation. c) The upper photograph shows the manufactured torsional oscillator. The lower photograph shows the complete setup including both permanent magnets and the torsional oscillator placed in the experimental well which was filled with the sample liquids.

Figure 5 :

 5 Figure 5: Torques considered in the principle of moment equilibrium. A double arrow illustrates an acting torque.

Figure 6 :

 6 Figure6: Comparison of the solutions for the fluid forces F c and F p acting on an oscillating cylinder and an in-plane oscillating plate, respectively. For both cases, the same surface interacting with the liquid is considered. In case of the solution for the cylinder, Hankel functions are obtained which make an intuitive interpretation of the obtained equations difficult. For high ratios of the cylinder's radius r c and the penetration depth δ, the difference between both solutions is negligibly small.

A c c e p t e d M a n u s c r i p t

  large amount of variables, which are usually not exactly known, constants and variables are combined in single factors where possible to reduce the amount of variables which have to be determined in a parameter fit.

Figure 7 :

 7 Figure 7: Rotor cross section

Figure 8 :

 8 Figure 8: Measured amplitudes as well as evaluated resonance frequencies and quality factors for two liquid series. The one dimensional shear wave model as well as the generalized model are fit to the measured data. The model parameters are given on the right hand side.

2 .

 2 Despite their limited accuracy, the mean values for f r and Q given in Tab. 1 are used to calculate the values for viscosity and mass density for both liquid series. The evaluation of absolute and relative errors given in Tab. 1 shows that with the present setup, absolute and relative accuracies for viscosity in the range of |∆ η| min , |∆ η| max = [0.003, 0.484] mPa•s and |∆ η rel | min , |∆ η rel | max = [0.003, 0.493], respectively are obtained. For the calculated mass densities the obtained accuracies are |∆ ρ| min , |∆ ρ| max = [0.0045, 0.0058] g/cm 3 and |∆ ρ rel | min , |∆ ρ rel | max = [0.005, 0.065]. This evaluation substantiates the need for setup improvement, especially to obtain

Figure 9 :

 9 Figure 9: Sensitivities to viscosity and mass density:The sensitivities of f r and Q to η and ρ are not constant but dependent both on η and ρ. For this reason the sensitivities are depicted as bands, in the experimentally investigated range of viscosities and mass densities. In the plots for S fr,η and S Q,η the upper boundary of the bands are the values for ρ = 1.01 g/cm 3 , the lower boundary for ρ = 0.79 g/cm3 . In the plots for S fr,ρ and S Q,ρ the upper and lower boundaries of the bands represent the evaluated values for η = 0.21 mPa•s and η = 2.05 mPa•s, respectively.

A

  c c e p t e d M a n u s c r i p t m 0k /s 2 m ρk / m 3 •s 2 kg m ηρk / m 2 kg • s 2 1.57 • 10 -7 4.64 • 10 -12 5.29 • 10 -8 c 0k /s c ηk / m•s 2 kg c ηρk / m 2 kg • s 2 1.88 • 10 -6 1.93 • 10 -4 6.73 • 10 -8

Figure 10 :

 10 Figure 10: Cross sensitivity of the resonance frequency to temperature. The upper figure shows the magnitudes of recorded frequency responses, the lower figure evaluated resonance frequencies versus temperature.

9

 9 well as evaluated resonance frequencies and quality factors for two liquid series. The one dimensional shear wave model as well as the generalized model are fit to the measured data. The model parameters are given on the right hand side. . . . . . . . . 7 Sensitivities to viscosity and mass density: The sensitivities of f r and Q to η and ρ are not constant but dependent both on η and ρ. For this reason the sensitivities are depicted as bands, in the experimentally investigated range of viscosities and mass densities. In the plots for S f r ,η and S Q,η the upper boundary of the bands are the values for ρ = 1.01 g/cm 3 , the lower boundary for ρ = 0.79 g/cm 3 . In the plots for S f r ,ρ and S Q,ρ the upper and lower boundaries of the bands represent the evaluated values for η = 0.21 mPa•s and η = 2.05 mPa•s, respectively. 8 10 Cross sensitivity of the resonance frequency to temperature. The upper figure shows the magnitudes of recorded frequency responses, the lower figure evaluated resonance frequencies versus temperature. . . . . . . . . . . . . . . . . . . . 9

Table 1 :

 1 Upper part: Acetone-isopropanol solutions for viscosity measurements. Lower part: Solutions for mass density measurements. The plus-minus values are evaluated typical errors (single standard deviations). x are the calculated values for viscosity and mass density using Eqs. 19 and 20 and evaluated mean values for f r and Q. ∆ x = xx and ∆ x rel = ∆ x/x are absolute and relative deviations from the values for viscosity and mass density, respectively.

	3 )	∆ ρ rel

Table 2 :

 2 Fitted model parameters using Eqs.19 and 20.

Acknowledgment

This work has been partially supported by the Linz Center of Mechatronics (LCM) in the framework of the Austrian COMET-K2 program. We also want to thank Bernhard Mayrhofer and Johann Katzenmayer for their help and excellent assistance.


E Decencière

E Tancrède-Bohin

P Dokládal

S Koudoro

A.-M Pena

T Baldeweck

Automatic 3D segmentation of multiphoton images: a key step for the quantification of human skin

Keywords: Multiphoton microscopy, 3D image processing, segmentation, graph cuts, watershed, quantification, human skin

Background/purpose: Multiphoton microscopy has emerged in the past decade as a useful non invasive imaging technique for in vivo human skin characterization. However, it has not been used until now in routine clinical practice, mainly because of the lack of specific image processing tools that would allow the investigator to extract pertinent quantitative three-dimensional (3D) information from the different skin components.

Methods: We propose a 3D automatic segmentation method of multiphoton images which is a key step for epidermis and dermis quantification. This method, based on the morphological watershed and graph cuts algorithms, takes into account the real shape of the skin surface and of the dermal-epidermal junction and allows separating in 3D the epidermis and the superficial dermis.

Results: The automatic segmentation method and the associated quantitative measurements have been developed and validated on a clinical database designed for aging characterization. The segmentation achieves its goals for epidermis-dermis separation and allows quantitative measurements inside the different skin compartments with sufficient relevance.

1

Conclusions: This study shows that multiphoton microscopy associated with specific image processing tools provides access to new quantitative measurements on the various skin components. The proposed 3D automatic segmentation method will contribute to build a powerful tool for characterizing human skin condition. To our knowledge, this is the first 3D approach to the segmentation and quantification of these original images.

In the past decade, multiphoton microscopy has emerged as a useful non invasive technique for in vivo skin imaging, as it allows 3D data (x-yz) on the epidermis and superficial dermis to be acquired with sub-cellular resolution. Multiphoton microscopy provides complement modalities, such as two-photon excited fluorescence (2PEF) and second harmonic generation (SHG), enabling non invasive imaging of skin to be performed through combined endogenous sources of contrast [START_REF] Baldeweck | In vivo multiphoton microscopy associated to 3D image processing for human skin characterization[END_REF][START_REF] Madani | In vivo multiphoton imaging of human skin: assessment of topical corticosteroid-induced epidermis atrophy and depigmentation[END_REF]. The specificity of these signals enables the various morphological structures of the skin to be identified. In epidermis, endogenous fluorophores (keratin, NAD(P)H, flavins, melanin, ...) contribute to the 2PEF signal. In dermis, two specific signals are generated: on the one hand, the 2PEF signal is produced by elastin, yielding images of the elastic fiber network. On the other hand the SHG signal originates from fibrillar collagens, allowing collagen fiber networks to be visualized. Therefore, this microscopy technique provides a powerful tool for characterizing superficial parts of skin and their follow-up under various conditions.

Until now, multiphoton microscopy has not been used in routine clinical practice, mainly because of the lack of specific image processing tools that would allow the investigator to extract pertinent quantitative 3Dinformation from the different skin components. Indeed, a very few image processing methods have been developed for characterizing multiphoton images of human skin in vivo. In dermis, the proposed methods are mainly based on the computation of mean intensity of 2PEF and SHG signals, or of the density of pixels occupied by these signals on conveniently chosen 2D slices within dermis at a fixed depth from the skin surface [START_REF] Koehler | In vivo assessment of human skin aging by multiphoton laser scanning tomography[END_REF][START_REF] Sugata | Evaluation of photoaging in facial skin by multiphoton laser scanning microscopy[END_REF]. More recently, parameters linked to the morphology of keratinocytes in epidermis have been evaluated in 2D on several slices [START_REF] Koehler | Keratinocyte morphology of human skin evaluated by in vivo multiphoton laser tomography[END_REF].

Assessing quantitative parameters in the 3D volume corresponding either to epidermis or to dermis requires to segment the multiphoton images in order to separate these compartments. Such a task is all but trivial, especially if the dermal-epidermal junction (DEJ) is neither flat nor parallel to the image plane. Up to our knowledge, no 3D approach to this issue has been previously proposed.

The goal of this work is to provide an automatic 3D segmentation method of multiphoton images of human skin, which is the key step for 3D quantification. This automatic 3D segmentation aims at separating the three main compartments appearing in the images: the coupling medium between skin and the cover glass (i.e. water), epidermis and dermis. Such a segmentation allows us to obtain a geometrical characterization of these compartments, as well as their interfaces, namely the skin surface and the dermal-epidermal junction. In addition, it opens the possibility to characterize the cells or the elastic and collagen fibers inside each 3D skin compartment in terms of quantity, size or organization.

Multiphoton skin images 1.Clinical study

The study involved 15 young (18-25 years) and 15 aged (70-75 years) human female volunteers. Before enrolment, each volunteer received information about the investigation procedure, was fully informed about the nature of the study and signed an informed consent. The experimental protocol was approved by the Saint-Louis Hospital (Paris) ethics committee and we followed the Declaration of Helsinki protocols. Multiphoton images were taken from the ventral side of the forearm.

Multiphoton setup -DermaInspect R device

Multiphoton imaging was performed with the DermaInspect R (JenLab GmbH, Jena, Germany), a CE-marked medical imaging system. SHG and 2PEF signals were simultaneously excited by a femtosecond Titanium-Sapphire laser adjusted to 760 nm as previously described [START_REF] Madani | In vivo multiphoton imaging of human skin: assessment of topical corticosteroid-induced epidermis atrophy and depigmentation[END_REF].

The acquisition protocol entailed a trade-off between the microscope characteristics (wavelength, objective), the optical properties of the living tissue, the size of the structures to be characterized, and practical considerations (mainly acquisition time due to in vivo imaging). As a result, the optical resolution was very anisotropic. The sampling step being chosen in accordance with this resolution, the images also were highly anisotropic (the sampling step in the z direction is 9 times larger than the sampling steps in the x and y directions).

In practice, each acquisition produces two 3D images of 511 × 511 × 70 voxels, with a sampling step of 0.255 µm in the plane orthogonal to the optical axis ((x, y) plane), and of 2.346 µm in the third direction (z axis). The first image corresponds to the 2PEF channel, and the second to the SHG channel. On each voxel, the grey level indicates the number of photons counted by the detectors.

Image description

The information contained in the 3D images is shown in Fig. 1. Top of the image (i.e. low z values) corresponds to the coupling water medium. There, respective values of the 2PEF and SHG channels are very low (< 4 photons), i.e. a noise range. Moving into the skin (increasing z) starts with the stratum corneum (Fig. 1a), composed of dead cells deprived of nuclei (corneocytes), and further the living epidermis (Fig. 1b-d), mainly containing keratinocytes. The lowest layer of epidermis is the basal layer (Fig. 1d), where keratinocytes divide before migrating upwards to the surface. This layer, also harboring melanocytes, typically shows the highest melanin density within keratinocytes. Logically these cells tend to generate a 2PEF signal stronger than the one produced by non melanin-enriched keratinocytes. Finally, going deeper through the DEJ interface, dermis is reached (Fig. 1e-f), where the 2PEF and SHG signals mainly originate from the elastic and collagen fibers.

Automatic segmentation method

The aim of the segmentation method is to classify each image voxel in one of the following compartments: coupling medium (CM), epidermis (ED) or dermis (D). The result will be a partition of the image domain into these three regions. Moreover, given that a typical clinical study based on these acquisitions may contain several hundred images, the segmentation method should not be time-consuming.

From a methodological image processing viewpoint, some characteristics of the multiphoton images have to be taken into account in the 3D segmentation method:

• The signal tends to be more intense near the central vertical axis of the image (i.e. the optical axis), and lower near the vertical borders.

• The signal intensity decreases with the imaging depth due to the attenuation and diffusion of excitation and emission light by the skin. • The vertical resolution and sampling steps are large with respect to some structures: for example, at the top of the living epidermis, cells are horizontal and flat, their thickness being possibly smaller than the vertical resolution.

• The sampling steps in the x-y-z directions are not isotropic: large sampling step in the z direction with respect to the sampling steps in the x and y directions.

• As the SHG signal is created by collagen fibers, one can use this information to identify the dermis. However, it is noteworthy that SHG signal alone does not allow the dermis to be delimited: at the DEJ level, dermal regions can be found with elastic fibers only and no collagen fibers, i.e. an absence of SHG signal.

All these features obviously make segmentation difficult. Based on the previous observations, we make the following assumptions the segmentation method will be based upon:

• The coupling medium (CM region) should be on the top of the image, and include the whole first slice. In this region, both 2PEF and SHG signals are low.

• The frontier between CM and ED regions, which should correspond to the skin surface, is characterized by a high gradient in the 2PEF signal. It tends to be horizontal.

• The upper part of epidemis (ED region), which should correspond to the stratum corneum, is indicated by high intensity 2PEF signal values, mainly generated by keratin in the corneocytes. In our experimental conditions, the mean 2PEF signal decreases when penetrating into the living epidermis, then increases again when coming near the DEJ. Indeed, the bottom layer of the epidermis, i.e. the basal layer, is often characterized by a high intensity 2PEF signal, at least near the optical axis.

• The upper part of dermis (D region) should theoretically show a significant SHG signal. However, as said above, collagen content may be highly variable among dermal areas near DEJ. On the other hand, 2PEF signal can be high in dermis when solar elastosis is present. Moreover, the junction between D and ED regions, corresponding to the DEJ, should be relatively smooth.

• Finally CM region should not be in contact with D region. The three regions can be seen as three layers, the one corresponding to the epidermis separating those corresponding to the coupling medium and the dermis.

The segmentation method will have to enforce these hypotheses.

Segmentation method overview

The segmentation strategy was designed to overcome the above mentioned difficulties.

A marker-based approach [START_REF] Meyer | Morphological segmentation[END_REF] has been chosen. Markers are sets of pixels considered as specific to a given compartment. They are afterwards propagated to the rest of the image domain, in order to have a complete partition. A marker-based approach is appropriate for our case, as the a priori information we have allows us to build them with a good level of confidence.

The propagation is done in two steps. Firstly, a fast watershed algorithm produces a low-level segmentation, made of small, homogeneous regions, also called super-voxels. Secondly, the super-voxels are merged, taking into account the markers, thanks to a graph cut algorithm.

The general strategy is summed up by the flow chart shown in Fig. 2. All processing is done in 3D, except otherwise specified.

Markers construction

The markers construction is based on thresholds of 2PEF and SHG signals.

The thresholds used to compute the coupling medium and epidermis markers are obtained by dividing the maximal value of the 2PEF histogram by empirically chosen constants. Similarly, the threshold used to compute the dermis marker is obtained by dividing the maximal value of the SHG histogram by a constant. The threshold is applied to the corresponding channel (low 2PEF values for the coupling medium; high 2PEF values for epidermis; high SHG values for dermis).

Conflicts (i.e. a non empty intersection) between markers of the coupling medium and epidermis do not occur, since the corresponding thresholds are mutually exclusive. Conflicts between epidermis and dermis markers may occur, since high 2PEF values are generated both in epidermis (melanin,..) and in dermis (elastin). In such a case, priority is given to the dermis marker, since the SHG signal theoretically only originates from dermis.

The connected components of the resulting images touching a reference plane are kept: the plane z = 0 for the coupling medium; z = 15 for the epidermis; and the plane with the highest mean SHG value for the dermis.

Markers propagation

Several methods have been tested to propagate the markers. The watershed algorithm based on markers appears very efficient, albeit segmentation errors frequently appear near the image vertical borders, given the lack of a significant gradient in those regions [START_REF] Beucher | Use of watersheds in contour detection[END_REF]. In fact, we gave privilege to a method which uses gradient information when available, and otherwise minimizes the frontier area. Energy minimization methods offer a good approach to the problem. Graph cuts belong to this category [START_REF] Boykov | Fast approximate energy minimization via graph cuts[END_REF]. A graph node is associated with each image voxel. When two voxels are neighboring, the corresponding nodes are linked with a graph edge. Two source nodes, traditionally called s and t, are added to the graph, and linked to the nodes which correspond to marked voxels. A weight is associated with each edge. The resulting weighted graph is finally partitioned into two sub-graphs, the first of which contains s, and the second t. This is done thanks to a flow minimization algorithm. The method is such that the sum of the weights of the removed edges is minimal.

Graph cuts have been applied to large images, with three dimensions or more [START_REF] Boykov | Graph cuts and efficient N-D segmentation[END_REF]. However, their theoretical complexity and memory requirements remain a problem in the present case, since several hundred 3D images have to be processed in a few hours. Therefore, instead of associating one node with each image voxel, we first compute a low level segmentation using the watershed approach, associate a node with each region of this low level segmentation, and add an edge between nodes such that the corresponding regions are neighbours. This approach was first proposed by Li et al. for 2D images [START_REF] Li | Lazy snapping[END_REF], and further developed and applied to 3D, single-modality images by Stawiaski and Decencière [START_REF] Stawiaski | Region merging via graph cuts[END_REF]. Here we extend it to the multi-modality case. This approach affords the possibility of computing edge weights which take more information into account.

This kind of graph cut can only segment a graph into two sub-graphs, despite the goal of segmenting the 3D image into three regions. Given such a layered model, two steps were carried out: the frontier between the coupling medium and the epidermis was computed first, starting from their respective markers and second, analogously, the frontier between epidermis and dermis was computed.

In practice, starting from the morphological gradient of the supremum between the 2PEF and SHG channels, the low level segmentation is computed in 2D, using a 4-connectivity, owing to the high anisotropy in the sampling steps (9 times larger in the z direction). The objective of the graph cuts can therefore be seen as the fusion of the regions of the low level segmentation. The last choice to be made is that of the weights of the graph edges, which is an essential aspect of the design of graph cuts [START_REF] Grady | Weights and topology: A study of the effects of graph construction on 3D image segmentation[END_REF]. Several functions were tested. Let us consider two neighboring regions A and B of the low level segmentation. N (A, B) denotes all the pairs of voxels (p, q) such that p and q belong to A and B respectively, and p and q are neighbors (a 3D 6-connectivity is considered). The weight w associated with the edge e linking the corresponding nodes is:

w(e) = (p,q)∈N (A,B) a(p, q) 1 + L 1 (I(p), I(q))
,

where:

• I(p) is the image value of voxel p, i.e. a couple (I 2P EF (p), I SHG (p)) corresponding to the two channels measured on that point,

• L 1 (I(p), I(q)) is the L 1 distance between the vector values I(p) and I(q), and

• a(p, q) is the area between voxels p and q.

Segmentation result

The final result of the segmentation is a 3D image containing three gray levels, corresponding to the coupling medium (dark gray in the illustrations), epidermis (medium gray) and dermis (light gray), respectively. As an example, Fig. 3 illustrates the result of the application of this 3D automatic segmentation method to two different cases: a young and an old volunteers.

The decrease in epidermis thickness and flattening of the dermal-epidermal junction with age are clearly evidenced. Fig. 4 shows the corresponding surface rendering of DEJ, whose flattening is clearly visible on these 3D reconstructions. These qualitative observations based on the segmentation could also be made on the raw data, however the segmentation opens the possibility of computing quantitative measurements, on the segmentation frontiers, and on each layer.

Quantification based on the 3D segmentation

Our automatic 3D segmentation of multiphoton skin images opens the possibility to characterize inside each 3D skin compartment the cells or the elastic and collagen fibers, in terms of quantity, size or organization.

Several measurements, based on this segmentation, have been developed and tested. Here we describe those that have been used to validate the automatic 3D segmentation method (see section 4) on the one hand, and to characterize aging (section 5), on the other hand. They are classified into three categories: morphological, density and organization measurements.

Morphological measurements

Morphological measurements are directly derived from the segmentation, without further need of the original signals. They characterize the regions and frontiers shapes. Here, we only consider two of them.

The epidermis thickness, T (ED), is the mean thickness, in µm, of the epidermis. The thickness of the epidermal compartment is measured along the z axis at each (x, y) position, and the mean value is computed.

The normalized area of DEJ is a generalization of the index of interdigitation described in two dimensions by Timár et al. on histological sections of human skin [START_REF] Timár | Interdigitation index -a parameter for differentiating between young and older skin specimens[END_REF]. Let A(DEJ) be the area of this surface, and A(H) the area of its projection on a horizontal plane. Then the normalized area is The more the DEJ is undulated, the larger the parameter value.

Density measurements

Sign 2P EF (D) and Sign SHG (D) give the proportion of significant pixels in the dermis, for the 2PEF and SHG images, respectively. A pixel is considered as significant if its value is higher than a given threshold, here equal to 5 photons, for both channels. This value is such that, given the noise present in the images, a signal value higher than this threshold almost surely corresponds to a pixel where the signal is not only the result of noise. These measurements are used to characterize the density of elastic and collagen fibers.

We also estimate in 3D a quantitative parameter characterizing the dermis, the SHG-to-AF aging index of dermis (SAAID), which has been calculated in 2D in the literature [START_REF] Koehler | In vivo assessment of human skin aging by multiphoton laser scanning tomography[END_REF][START_REF] Sugata | Evaluation of photoaging in facial skin by multiphoton laser scanning microscopy[END_REF]. 3D-SAAID of the dermis is a straightforward generalization of 2D-SAAID:

3D -SAAID(D) = Sign SHG (D) -Sign 2P EF (D) Sign SHG (D) + Sign 2P EF (D) .
3D-SAAID(D) is between -1 and 1. Negative values mean that there are more significant pixels with 2PEF than with SHG signals, and a value close to zero means that both proportions are similar. However, this value says nothing about their relative positions: do they coincide in the same pixels or not? In the following section, we propose a new measurement for this purpose.

New multimodal organization measurement

Multimodal organization measurements aim at quantifying the relative position of 2PEF and SHG signals. In order to measure the coincidence of 2PEF and SHG, we propose to measure the normalized overlapping N OL(D) of both signals in the dermis:

N OL(D) = D 2P EF n SHG n D 2P EF n SHG n
, where and are the infimum and supremum operators, respectively, and 2P EF n and SHG n correspond to the 2PEF and SHG signals, normalized on the dermis so that their integrals are equal, i.e.: D 2P EF n = D SHG n . N OL(D) is between 0 and 1; a large value means that the signals are spatially correlated; a low value indicates that they tend to be mutually exclusive. To our knowledge, this is the first time that this measurement is applied to multiphoton human skin images.

Validation of the segmentation method

Validating an automatic 3D segmentation method on complex biological data is a difficult task. The most straightforward approach is to ask at least two experts to produce a reference segmentation on a representative number of images.

Reference segmentation

Providing a reference segmentation on a representative number of images is an extremely tedious task. Therefore, we choose to build for each 3D multiphoton image of the data base a reference segmentation obtained by manually correcting the automatic segmentation on the whole data base. For that purpose, two experts inspected all the automatic segmentations, and manually corrected those which were not totally satisfactory. The resulting segmentations are hereafter called, for simplicity, the manual segmentations. Finally, measurements were recomputed using the manual segmentations.

Performance evaluation

One has to keep in mind that segmentation is only an intermediate step. The final goal of our clinical studies is to obtain an automatic segmentation, not necessarily pixel-to-pixel correlated to the manual one, but leading to a relevant estimation of the different quantitative measurements. Therefore, the issue is: are the measurements obtained thanks to these automatic segmentations accurate enough? For answering this question, the measurements described in the previous section were applied to the result of the automatic segmentation. If the proposed segmentation method is good enough for the current application, then the data obtained on the automatic and manual segmentations should be correlated, so that the same conclusions can be drawn from both sets of measurements. Fig. 5 shows the correlation graphs for the six described measurements. It appears that the correlation does seem good. Indeed correlation coefficients are systematically above 0.84. Correlation coefficients of morphological measurements, i.e. epidermis thickness (T (ED)) and normalized area of DEJ (N A(DEJ)) are equal to 0.88 and 0.84, respectively; they are smaller than those of the four other measurements within the dermis, which are higher than 0.93. As expected, the quality of segmentation has more impact on the measurements linked to the shape of the frontiers delimiting the epidermis (i.e. T (ED) and N A(DEJ)), than on those linked to the content of the skin compartments. Accordingly, one should be careful when using this automatic segmentation method in studies where a slight shape alteration of DEJ is to be characterized. In this case, a manual correction of the segmentation will be necessary.

Concerning the linear regression (as displayed on Fig. 5), the equation of the linear regression for each measurement leads to a slope coefficient close to 1, and to an ordinate at origin close to 0, which is a sign of non-bias.

In conclusion, we can state that our automatic segmentation achieves its goals and allows quantitative measurements inside the different skin compartments with sufficient relevance. Moreover, the current version of the segmentation method, which has not undergone any specific computational optimization, takes an average of one minute to process a 3D multiphoton image on a classical personal computer (Intel core i7, 8 GB memory).

Application to aging characterization

The clinical study (see section 1) has been designed to evidence age-related (young/old) differences. The aim of the study was to evaluate the relevance and interest of the 3D measurements obtained with the automatic segmentation method. Resulting data is given in Fig. 6.

Morphological measurements

For the mean thickness of epidermis, T (ED), a decrease with age is noticed: mean value in the young group is 66 µm, whilst it decreases to 59 µm in the old group.

Similarly, the normalized area of DEJ, N A(DEJ), which is related to the shape of DEJ, decreases with age (2.40 in the young group; 2.05 in the old group).

These age-dependent morphological changes are expected: with age, the epidermis becomes thinner and the rete ridges of the dermal-epidermal junction flatten out. It should be noted that, given the relatively small field of view, it was not evident that the measurements would reveal such morphological changes. Our 3D results are in agreement with the histology and other imaging techniques (optical coherence tomography, ultrasound imaging) performed on the same volunteers (data not shown) as well as with anual en x ? published 2D data obtained with histology [START_REF] Lavker | Aged skin: a study by light, transmission electron, and scanning electron microscopy[END_REF] or with other in vivo imaging techniques (optical coherence tomography, confocal microscopy) [START_REF] Sauermann | Age related changes of human skin investigated with histometric measurements by confocal laser scanning microscopy in vivo[END_REF][START_REF] Querleux | Skin from various ethnic origins and aging: an in vivo cross-sectional multimodality imaging study[END_REF].

Density measurements As regards Sign 2P EF value, which is linked to the density of elastic fibers, an increase with aging is found. The results are in agreement with the histology performed on the same volunteers (data not shown). Sign SHG value, which is linked to the density of collagen fibers, shows no difference between old and young volunteers.

We have also found that 3D-SAAID highly decreases with aging, as previously shown in 2D [START_REF] Lin | Evaluating cutaneous photoaging by use of multiphoton fluorescence and second-harmonic generation microscopy[END_REF][START_REF] Koehler | In vivo assessment of human skin aging by multiphoton laser scanning tomography[END_REF][START_REF] Sugata | Evaluation of photoaging in facial skin by multiphoton laser scanning microscopy[END_REF][START_REF] Koehler | Intrinsic, solar and sunbed-induced skin aging measured in vivo by multiphoton laser tomography and biophysical methods[END_REF]. This difference is directly linked to the relative increase in significant 2PEF pixels versus significant SHG pixels.

New organization measurement Lastly, the normalized overlapping N OL(D) in the dermis tends to increase with age. This means that the correlation between 2PEF and SHG signals increases with age, which might correspond to a more disordered structure. This new parameter might be a new quantitative descriptor of the modifications of the superficial dermis organization during aging.

All the measurements have been made on the whole 3D region, either corresponding to epidermis or dermis in the 3D image. In our opinion, this is a clear advantage with respect to what has been previously reported in the literature for the dermis, i.e. taking empirically a 2D slice which is completely included within the dermis [START_REF] Koehler | In vivo assessment of human skin aging by multiphoton laser scanning tomography[END_REF][START_REF] Kaatz | Depth-resolved measurement of the dermal matrix composition by multiphoton laser tomography[END_REF][START_REF] Sugata | Evaluation of photoaging in facial skin by multiphoton laser scanning microscopy[END_REF]. Not only the present approach avoids the selection of a 2D plane, but the measurements are based on a larger, more representative, volume of tissue.

In addition, all the data from automatic segmentation is consistent with data based on the manually corrected segmentation, which confirms that our automatic segmentation method is suitable for reliable quantitative characterization of human skin from multiphoton images.

Conclusion and perspectives

A 3D automatic segmentation method for analyzing in vivo multiphoton images of human skin has been presented. To our knowledge, this is the first 3D approach to the quantification of these images. Our method takes into account the real shape of the dermal-epidermal junction and is able to deal both with anisotropic sampling steps and spatially varying signal to noise ratio. Existing and new measurements have been applied in 3D inside each segmentation region. The automatic segmentation method and the associated quantitative measurements have been developed and validated on a clinical database designed for aging characterization. It is shown that the segmentation achieves its goals, and that the obtained measures are in agreement with previously reported results.

Segmentation is a key step prior to 3D quantification of the structures within each compartment, like epidermal cells and dermal fibers. This work has prompted several ongoing studies, such as the correction of the multiphoton signal attenuation, the development of new parameters for epidermis and dermis characterization, and the evaluation of the effects of different cosmetic and medical products on the skin using this method.

This study shows that multiphoton microscopy associated with specific image processing tools provides access to new relevant quantitative measurements on the various skin components. The development of this multidisciplinary approach associating multiphoton microscopy and image processing tools will be essential for the future deployment of this imaging technique in routine clinical practice.

Figure 1 :

 1 Figure 1: In vivo multiphoton images of normal human skin. 2PEF signal (cyan hot color) reveals the endogenous fluorophores distribution inside (a) corneocytes in the stratum corneum, keratinocytes in the upper (b) and lower (c) living epidermis, (d) keratinocytes in the stratum basale. In the dermis the 2PEF signal reveals elastic fibers whereas the SHG signal (red color) reveals the collagen fiber distribution: near DEJ fibers are clearly visible (e), whereas they become fuzzier with depth (f). Acquisition conditions: 40×, 1.3-NA oil immersion objective; excitation power exponentially increased from 12 mW at the surface of the skin up to 47 mW at a depth of 75 µm and kept constant for imaging depths greater than 75 µm; acquisition time: 7.4 seconds per image of 511 × 511 pixels and 8.5 minutes per stack of 130 × 130 µm 2 × 70 images acquired every 2.346 µm.

Figure 2 :

 2 Figure 2: Flow chart of the 3D automatic segmentation method.

Figure 3 :Figure 4 :

 34 Figure 3: 3D results of the automatic segmentation applied to multiphoton images from a young (top) and an old (bottom) volunteers: (a) and (a') raw 2PEF/SHG data; (b) and (b') 3D automatic segmentation (coupling medium -dark gray, epidermis -medium gray and dermis -light gray); (c) and (c') overlay of 3D automatic segmentation on raw data.

Figure 5 :

 5 Figure 5: Correlation graphs for the 6 measurements: T (ED), N A(DEJ), Sign 2P EF (D), Sign SHG (D), 3D-SAAID(D) and N OL(D). Horizontal axis: values obtained with automatic segmentation; vertical axis: values obtained with manual segmentation. Each blue point represents the result obtained for one 3D multiphoton image; the linear regression fit is shown with a red line.

Figure 6 :

 6 Figure 6: Quantification results obtained on the aging database: change with age in the mean thickness of epidermis (a), the normalized area of DEJ (b), the density of elastic (c) and collagen (d) fibers, the SHG-to-AF aging index of dermis (e) and the normalized overlapping of 2PEF and SHG signals in the dermis (f). Error bars correspond to 95% confidence intervals of the mean.


Closed-forms of Kirchhoff elastic rods shape and sensitivity in the planar case

Olivier Roussel 1 , Marc Renaud 1 and Michel Taïx 1Å

bstract

In this report we give closed-forms of Kirchhoff 3-D elastic rods curvature in terms of elliptic functions and, by treating planar rods as a special case, we show we can also obtain closed-forms of planar rods shape, sensitivity and total elastic energy.

General case of 3-D rods

Consider an inextensible, non-shearable and unit length linearly elastic rod. The shape of the rod traces a curve that we will describe by the mapping q : r0, 1s Ñ SEp3q. The position along the rod is parametrized by t P r0, 1s and we will name "base" and "tip" of the rod its extremity at t " 0 and t " 1 respectively. Let the mappings u 1 ptq, u 2 ptq, u 3 ptq such that u i : r0, 1s Ñ R be axial and bending rod strains respectively, and c 1 , c 2 , c 3 be the constants that reflect its elasticity properties. As in [START_REF] Jurdjevic | Integrable Hamiltonian systems on complex Lie groups[END_REF], we say the elastic rod is in static equilibrium in the sense of Kirchhoff if it locally minimizes the elastic energy defined by

E el " 1 2
ż 1 0 3 ÿ i"1 c i u 2 i dt. (1.1)
Without loss of generality, we will also assume that the base of the rod is held fixed at the origin, i.e. qp0q " e where e is the identity element of SEp3q. Under these assumptions, we will denote by B the set of positions that the other extremity of the rod qp1q can reach. As shown in [START_REF] Bretl | Quasi-static manipulation of a kirchhoff elastic rod based on a geometric analysis of equilibrium configurations[END_REF], the problem of static equilibrium of such rods can be formulated as an optimal control problem by minimize q,u 1 2

ż 1 0 3 ÿ i"1 c i u 2 i dt subject to 9 q " q ˜3 ÿ i"1 u i X i `X4 qp0q " e, qp1q " b (1.2)
for some b P B and where

X 1 " " 0 0 0 0 0 0 ´1 0 0 1 0 0 0 0 0 0  X 2 " " 0 0 1 0 0 0 0 0 ´1 0 0 0 0 0 0 0  X 3 " " 0 ´1 0 0 1 0 0 0 0 0 0 0 0 0 0 0  X 4 " " 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0  X 5 " " 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0  X 6 "
" 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0  is a basis for sep3q, the Lie algebra of SEp3q. Note that when solving this optimal control problem, the rod tip position b is not an input.

In these conditions, the Maximum Principle states that solutions to this optimal control problem are the projections of extremal curves defined on the cotangent bundle T ˚SEp3q onto SEp3q. Thanks to the Lie Group structure of SEp3q, the Hamiltonian can be reduced on the dual of the Lie algebra sep3q ˚and the corresponding (time-varying) Hamiltonian vector fields µ : r0, 1s Ñ sep3q ˚can be expressed by

$ ' ' ' ' ' ' ' & ' ' ' ' ' ' ' % 9 µ 1 " µ3µ2 c3 ´µ2µ3 c2 9 µ 2 " µ 6 `µ1µ3 c1 ´µ1µ3 c3 9 µ 3 " ´µ5 `µ1µ2 c2 ´µ1µ2 c1 9 µ 4 " µ3µ5 c3 ´µ2µ6 c2 9 µ 5 " µ1µ6 c1 ´µ3µ4 c3 9 µ 6 " µ2µ4 c2 ´µ1µ5 c1 (1.3)
where vector fields µ are related to controls u i by u i " c ´1 i µ i for i P t1, 2, 3u. Let A be the set homeomorphic to R 6 and a P A such that a i fi µ i p0q, i P t1, . . . , 6u. It has been shown in [START_REF] Bretl | Quasi-static manipulation of a kirchhoff elastic rod based on a geometric analysis of equilibrium configurations[END_REF] that coordinates in A offer a global parameterization to the set of static equilibrium configuration for the rod. In other words, we can describe configurations of quasi-static 3-D elastic rods using the 6-dimensional configuration space A.

Assuming isotropy and normalized elasticity constants such that c i " 1 for i P t1, 2, 3u, we have from (1.3) 9 µ 1 " 0. Then µ 1 is a constant of motion with µ 1 " a 1 and

$ ' ' ' ' & ' ' ' ' % 9 µ 2 " µ 6 9 µ 3 " ´µ5 9 µ 4 " µ 3 µ 5 ´µ2 µ 6 9 µ 5 " a 1 µ 6 ´µ3 µ 4 9 µ 6 " µ 2 µ 4 ´a1 µ 5 (1.4)
The signed curvature κ and the torsion τ of the curve can be expressed in terms of µ by

κ 2 " µ 2 2 `µ2 3 τ " µ 1 ´µ2 µ 5 `µ3 µ 6 µ 2 2 `µ2 3
and, as mentioned in [START_REF] Bretl | Quasi-static manipulation of a kirchhoff elastic rod based on a geometric analysis of equilibrium configurations[END_REF], the differential system (1.4) is equivalent to 2: κ `κ3 ´2κpτ ´λ1 q 2 " λ 2 κ (1.5a)

κ 2 pτ ´λ1 q " λ 3 (1.5b)
where the constants of integration are given by

$ & % λ 1 fi a1 2 λ 2 fi a 2 2 `a2 3 `2a 4 ´a2 1 2 λ 3 fi a1 2 pa 2 2
`a2 3 q ´pa 2 a 5 `a3 a 6 q.

Substituting (1.5b) into (1.5a) and integrating, we obtain

9 κ 2 `1 4 κ 4 `λ2 3 κ ´2 ´λ2 2 κ 2 " λ 4 (1.6)
where the constant of integration λ 4 is given by

λ 4 fi a 2 5 `a2 6 ´1 4 pa 2 2 `a2 3 q 2 `1 2 pa 2 2 `a2 3 qpa 2 1 ´2a 4 q ´a1 pa 2 a 5 `a3 a 6 q
By making the change of variable υ " κ 2 , (1.6) transforms to

9 υ 2 `υ3 ´2λ 2 υ 2 ´4λ 4 υ `4λ 2 3 " 0 (1.7)
Figure 1: Plot of the cubic polynomial P pυq with respect to the squared curvature υ. Hatched regions correspond to impossible values illustrating that the only valid range for υ is given by α 2 and α 3 , the zeros of P pυq.

As already stated in [START_REF] Langer | The total squared curvature of closed curves[END_REF], this equation is in the form 9 υ 2 " P pυq with P is the cubic polynomial

P pυq fi ´υ3 `2λ 2 υ 2 `4λ 4 υ ´4λ 2 3 . (1.8)
Let ´α1 , α 2 , α 3 be the zeros of the polynomial P pυq such that ´α1 ď 0 ď α 2 ď α 3 .

(1.9)

. As P p˘8q " ¯8 and P p0q " ´4λ 2 3 ď 0, P pυq is in the form illustrated in figure 1. Also, we have υ ě 0 and P pυq ě 0 as they are both squares, so υ P rα 2 , α 3 s. The polynomial P pυq can be rewritten for its zeros by P pυq " ´pυ `α1 qpυ ´α2 qpυ ´α3 q.

We can express the polynomial zeros ´α1 , α 2 , α 3 from the constants of integrations λ i by

α 1 ´α2 ´α3 " ´2λ 2 α 1 α 2 `α1 α 3 ´α2 α 3 " 4λ 4 α 1 α 2 α 3 " 4λ 2 3 .
(1.10)

The squared curvature υ can be expressed in terms of elliptic functions by

υptq " α 3 ´1 ´n sn 2 `rt `ϕ|m ˘¯(1.11)
the parameter m, the characteristic n and r can be expressed from the polynomial zeros by

m " α 3 ´α2 α 3 `α1 n " α 3 ´α2 α 3 r " 1 2 ? α 3 `α1 (1.12) Given fi d 1 n ˆ1 ´a2 2 `a2 3 α 3 ˙(1.13)
the phase ϕ can be retrieved from a 2 2 `a2 3 " α 3 `1 ´n sn 2 pϕ|mq ˘, and is given by ϕ " sgnpa 3 a 5 ´a2 a 6 q arcsnp |mq (1.14

)
where arcsn is the inverse of the Jacobi elliptic function sn.

Note that from (1.9), we have 0 ď m ď n ď 1.

As outlined in [START_REF] Jurdjevic | Integrable Hamiltonian systems on complex Lie groups[END_REF], it has been shown the Hamiltonian vector fields in (1.4) is integrable and we have proved it can be expressed in the following form

$ ' ' ' ' ' & ' ' ' ' ' % µ 2 " κ sin ψ µ 3 " κ cos ψ µ 4 " 1 2 pλ 2 `a2 1 2 ´υq µ 5 " ´9 κ cos ψ `κ 9 ψ sin ψ µ 6 " 9 κ sin ψ `κ 9 ψ cos ψ.
(1.15)

where

ψptq " λ 1 t ´λ3 α 3 r ˆΠ´n , am `rt `ϕ|m ˘|m ¯´Π ´n, am `ϕ|m ˘|m ¯˙`ψ p0q
with Πpn, u|mq the elliptic integral of the third kind and ampu|mq is the Jacobi amplitude.

Planar case

Although neither the curve qptq nor the rod sensitivity Bqptq Ba can be explicitly expressed in the general 3-D case, we will show in this section that closed forms can be obtained in the planar case which can be treated as a particular case of the previously presented model.

Curvature and internal wrenches

Considering only planar curves qptq in the xy-plane with q " p0, 0, θ, x, y, 0q T , Hamiltonian vector fields defined in (1.3) simplify to

$ ' ' ' ' ' ' & ' ' ' ' ' ' % 9 µ 1 " 0 9 µ 2 " 0 9 µ 3 " ´µ5 9 µ 4 " µ 3 µ 5 9 µ 5 " ´µ3 µ 4 9 µ 6 " 0 (2.1)
Closed-forms of rod internal wrenches µptq defined in (1.15) reduce to

$ ' ' ' ' ' ' & ' ' ' ' ' ' % µ 1 " 0 µ 2 " 0 µ 3 " κ µ 4 " ´1 2 pκ 2 `λ2 q µ 5 " ´9 κ µ 6 " 0 (2.2)
And constants of integration defined in (1) simplify to

$ ' ' & ' ' % λ 1 " 0 λ 2 " a 2 3 `2a 4 λ 3 " 0 λ 4 " a 2 5 ´a2 3 p 1 4 a 2 3 `a4 q
We retrieve the same results as we would have obtain by applying the same problem formulation on the Lie Group SEp2q rather than SEp3q. Therefore, in the rest of this section we will restrict to solutions of (1.2) that are similar to trajectories on SEp2q, which are generated by the subset of initial conditions ta P A : pa 1 , a 2 , a 6 q " p0, 0, 0qu.

In the following equations, when referring to an elliptic function pq, we will simplify the notation pqpu|mq to pq u. Also, let us define

Γptq fi rt `ϕ
and the following constants of motion that be needed in the following developments

ε fi sgnpa 3 q Γ 0 fi Γp0q δ fi λ 2 2 `4λ 4
The expression of the phase ϕ given in (1.14) simplifies to ϕ " sgnpa 3 a 5 q arcsnp |mq where given in (1.13) reduces to

" d 1 n ˆ1 ´a2 3 α 3 ˙.

Expression of the curvature

In the planar case, as λ 3 " 0, the polynomial P pυq simplifies to P pυq " ´υ3 `2λ 2 υ 2 `4λ 4 υ so P pυq has one trivial zero at υ " 0.

From (1.9), we can distinguish three cases as outlined in [START_REF] Langer | The total squared curvature of closed curves[END_REF] and [START_REF] Singer | Lectures on elastic curves and rods[END_REF]:

• Case I: λ 4 ą 0
Using (1.10), we have that λ 1 pλ 2 `λ3 q ą λ 2 λ 3 .

This imposes the choice for the zeros to

$ & % α 1 " ´λ2 `?δ α 2 " 0 α 3 " λ 2 `?δ.
From (1.12) we get n " 1, so the squared curvature formula in (1.11) simplifies to υptq " α 3 `1 ´sn 2 Γptq "

α 3 cn 2 Γptq.
Then the signed curvature is given by κptq " ε ? α 3 cn Γptq.

(2.

3)

The curvature κptq oscillates between ? α 3 and ´?α 3 and the resulting curve qptq is called a "wavelike" elastica.

• Case II: λ 4 ă 0 Using (1.10), we have

λ 1 pλ 2 `λ3 q ă λ 2 λ 3 .
This imposes the choice for the zeros to

$ & % α 1 " 0 α 2 " λ 2 ´?δ α 3 " λ 2 `?δ.
From (1.12) we get n " m, so the squared curvature formula in (1.11) simplifies to υptq " α 3 `1 ´m sn 2 Γptq "

α 3 dn 2 Γptq.
Then the signed curvature is given by

κptq " ε ? α 3 dn Γptq (2.4)
The curvature κptq is non-vanishing and the resulting curve qptq is called a "orbit-like" elastica.

• Case III: λ 4 " 0 This borderline case implies the polynomial P pυq reduces to

P pυq " ´υ3 `2λ 2 υ 2
which has a double zero.

Using first equation of (1.10), only one choice is possible for the zeros α i :

" α 1 " α 2 " 0 α 3 " |2λ 2 | .
which leads to the signed curvature

κptq " ε ? α 3 sech Γptq (2.5)
This corresponds to the borderline case where the curvature is non-periodic.

Reduction to a unique formulation of the curvature

These cases can be reduced to a single formulation of the curvature by allowing the parameter m to be any positive or null real and applying the Jacobi's real transformation (see [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF] §16.11).

By relaxing the constraint on the zeros α i given in (1.9), and keeping only one fixed choice on the zeros that we will denote by α 1 such that

$ & % α 1 1 " 0 α 1 2 " λ 2 ´?δ α 1 3 " λ 2 `?δ.
In this form, α 1 3 is positive as

? δ ą |λ 2 | but α 1
2 can now be negative. Note that we still have α 1

3 ě α 1 2 . Using same forms as in (1.12), the elliptic parameter m 1 and r 1 by

m 1 " α 1 3 ´α1 2 α 1 3 r 1 " 1 2 a α 1
3
but as mentioned before, the new elliptic parameter m 1 is only constrained to in r0, 8q. Then, the signed curvature can be expressed by a unique expression by

κptq " ε a α 1 3 dn `r1 pt `φq ˇˇm 1 ˘(2.6)
When m 1 ą 1, the Jacobi's real transformation can be applied to reduce to a parameter m such that 0 ď m ď 1 and we retrieve the previously described cases.

Explicit formulation of rod total elastic energy

Recall from (1.1) the total elastic energy of the rod is given by

E el " 1 2 ż 1 0 u 3 ptq 2 dt " 1 2
ż 1 0 κptq 2 dt
Using the unique formulation of the curvature κptq given in (2.6) can be integrated to give an explicit formulation in terms of the elliptic integral of the second kind by

E el " α 1 3 2 E `r1 pt `φq ˇˇm 1 2.2

Integration of the curve qptq

From the differential system defined in (1.2), it follows that 9 θ " u 3 " κ 9

x " cos θ 9 y " cos θ.

Using (2.2), the integration of the curvature is given by

cos θptq " β 1 p0qβ 1 ptq `4β 2 p0qβ 2 ptq (2.7a) sin θptq " 2 ε pβ 1 p0qβ 2 ptq ´β2 p0qβ 1 ptqq (2.7b) xptq " β 1 p0q ż β 1 ptq `4β 2 p0q ż β 2 ptq (2.7c) yptq " 2 ε ˆβ1 p0q ż β 2 ptq ´β2 p0q ż β 1 ptq ˙. (2.7d)
The functions β 1 ptq and β 2 ptq can be explicitly given using Jacobi elliptic functions and the elliptic integral of second kind Epu|mq in the three cases previously described as follows

• Case I: λ 4 ą 0
Integrating the curvature in (2.3) (see [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF] §16.24) leads to θptq " 2 ε parccos pdn Γptqq ´arccos pdn Γ 0 qq Let Aptq fi arccos pdn Γptqq and Ap0q fi arccos pdn Γ 0 q, then cos Aptq " dn Γptq sin Aptq " ˘b1 ´dn

Explicit formulation of elastic rod sensitivity

In the 3-D case, the elastic rod sensitivity is given by the 6-dimensional Jacobian matrix where ˚represents indeterminate values.

As we can only obtain closed-forms of the rod shape and thus of rod sensitivity in this special, we will focus in this section on the 3-dimensional block J P pt, aq of Jpt, aq for i, j P t3, 4, 5u. Differentiating the general form of the curve qptq in (2.7) Regardless the three cases of curve elastica, we can derivate with respect to a the following forms:

  Γptqq ´arcsin psn Γ 0 qq Let Aptq fi arcsin psn Γptqq and Ap0q fi arcsin psn Γ 0 q, then ´2 dn 2 Γptq ˘`1 ´2 dn 2 Γ 0 ˘`4 cn Γptq sn Γptq cn Γ 0 sn Γ 0 ´2 sn 2 Γ 0 ˘cn Γptq sn Γptq ´`1 ´2 dn 2 Γptq ˘cn Γ 0 sn Γ 0 which is in the form (2.7) with β 1 and β 2 given byβ 1 ptq fi 1 ´2 sn 2 Γptq

	Given that θptq 2 " ε pAptq ´Ap0qq, we have Using half-angle formulas, we get
	cos cos θptq " cos 2 θptq 2 ´sin 2 θptq θptq " cos Aptq cos Ap0q `sin Aptq sin Ap0q 2 2 " cn Γptq cn Γ 0 `sn Γptq sn Γ 0 " `2 sech 2 Γptq ´1˘`2 sech 2 Γ 0 ´1˘`4 sech Γptq tanh Γptq sech Γ 0 tanh Γ 0
	sin θptq " 2 cos	sin θptq 2 sin θptq 2 θptq " sin Aptq cos Ap0q ´cos Aptq sin Ap0q 2
	" 2 ε ``2 sech 2 Γ 0	2 Γptq " ε psn Γptq cn Γ 0 ´cn Γ 0 sn Γptqq ´1˘s ech Γptq tanh Γptq ´`2 sech 2 Γptq	´1˘s ech Γ 0 tanh Γ 0	which
	" is in the form (2.7) with β 1 and β 2 given by ? m sn Γptq Using half-angle formulas, we get
	Given that θptq 2 " ε pAptq ´Ap0qq, we have cos θptq " cos 2 θptq 2 β 1 ptq fi 2 sech 2 Γptq ´1 ´sin 2 θptq 2 β 2 ptq fi sech Γptq tanh Γptq		(2.12a) (2.12b)
	cos sin Using half-angle formulas, we get θptq 2 " cos Aptq cos Ap0q `sin Aptq sin Ap0q " dn Γptq dn Γ 0 `m sn Γptq sn Γ 0 θptq 2 " sin Aptq cos Ap0q ´cos Aptq sin Ap0q " ε ? m psn Γptq dn Γ 0 ´sn Γ 0 dn Γptqq " `1 sin θptq " 2 cos θptq 2 which integrate to θptq ż sin 2 β 1 ptq " 2r ´1 ptanh Γptq ´tanh Γ 0 q ´t ż " 2ε ``1 (2.10a) (2.13a) β 2 ptq " ´r´1 psech Γ ptq ´sech Γ 0 q . (2.13b)
	cos θptq " cos 2 θptq 2 " `2 dn 2 Γ 0 and can be integrated to ´sin 2 θptq β 2 ptq fi sn Γptq cn Γptq 2 ´1˘`2 dn 2 Γptq ´1˘`4 m dn Γptq sn Γptq dn Γ 0 sn Γ 0 ż	(2.10b)
	sin θptq " 2 cos β 1 ptq " m θptq 2 " 2 ε ż ? m ``2 dn 2 Γ 0 θptq sin 2 β 2 ptq " ´r´1 pcn Γptq ´cn Γ 0 q ´1˘d n Γptq sn Γptq ´`2 dn 2 Γptq	´1˘d n Γ 0 sn Γ 0 (2.11b) which
	is in the form (2.7) with β 1 and β 2 given by		(2.11c)
			β 1 ptq fi 2 dn 2 Γptq ´1 β 2 ptq fi ? m sn Γptq dn Γptq		(2.8a) (2.8b)
	and can be integrated to
		ż
		β 1 ptq " 2r ´1 pEpam Γptqq ´Epam Γ 0 qq ´t		(2.9a)
		ż
		β 2 ptq " ´r´1 pcn Γptq ´cn Γ 0 q cos Aptq " ˘`1 `sinh 2 Γptq ˘´1 2		(2.9b)
			" sech Γptq		(2.9c)
	• Case II: λ 4 ă 0		sin Aptq " tanh Γptq
	Integrating the curvature in (2.4) (see [1] §16.24) leads to
	Given that θptq 2 " ε pAptq ´Ap0qq, we have cos θptq 2 " cos Aptq cos Ap0q `sin Aptq sin Ap0q θptq " 2 ε parcsin psn cos Aptq " ˘a1 ´sn 2 Γptq " sech Γptq sech Γ 0 `tanh Γptq tanh Γ 0
			" cn Γptq
			sin Aptq " sn Γptq

´1

`t pm ´2q `2r ´1 pE pam Γ ptqq ´E pam Γ 0 qq ˘(2.11a)

• Case III: λ 4 " 0 Integrating the curvature in (2.5) (see

[START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF]

§16.24) leads to θptq " 2 ε parctan psinh Γptqq ´arctan psinh Γ 0 qq Let Aptq fi arctan psinh Γptqq and Ap0q fi arctan psinh Γ 0 q, then sin θptq 2 " sin Aptq cos Ap0q ´cos Aptq sin Ap0q

" ε ptanh Γptq sech Γ 0 ´tanh Γ 0 sech Γptqq

Ṫhen, most of these forms simplify in the three previously introduced cases and we can give the explicit forms of the derivatives of functions β 1 ptq and β 2 ptq (and their respective integrals) with respect to a: