Laetitia Pernod

Antoine Ducoin

Hervé Le Sourne

Jean-François Sigrist

Coupled Numerical Simulation of an Aluminum and a Composite Hydrofoil in Steady and Unsteady Flows

Keywords: Composite hydrofoil, fluid-structure interaction, commercial software, numerical simulation, experimental validation

The objective of this paper is to investigate the fluid structure interaction of a 3D NACA0009 composite hydrofoil in steady and unsteady flows. This hydrofoil set up, which was originally developed experimentally by the Australian Maritime College (AMC), is clamped at the root and free at the tip section, and was designed to behave like a marine propeller blade (Zarruk et. al, 2014). 3D fluid structure computations are performed using a coupling between CFD (Computational Fluid Dynamic) and CSD (Computational Solid Dynamic) softwares, through tightly coupled algorithm. The hydrodynamic performances, structure natural frequencies and tip displacements are validated against experiments. Then the hydroelastic response is investigated. The results will help for the validation of a tightly coupled procedure, which will be used for the development and the optimization of a newly developed design of composite marine propeller.

INTRODUCTION

The development of composite marine propellers for surface ship and submarines has been the subject of increasing researches in the last decade. Within the main objective of CO2 emissions reduction, the key design clearly remains in the development of lighter structures, but also in the modification of the structural properties to increase the propeller performances and to reduce the possible vibrations caused by the flow. Academic studies demonstrated that using a passive control of composite propellers thanks to their structural properties can successfully lead to performance improvements and reduced vibrations [START_REF] Plucinski | Optimization of a Self-Twisting Composite Marine Propeller using a Genetic Algorithms[END_REF], (He et. al, 2012). As a consequence, there is now a growing interest for Navies and shipyards to develop full scale prototypes [START_REF] Yamatogi | Study on Cavitation Erosion of Composite Materials for Marine Propeller[END_REF].

It was shown that using composite materials instead of the usual metal alloys for marine propellers would be a suitable solution to face the challenge. The passive control of deformations typically allows to reduce cavitation [START_REF] Chen | Design, Fabrication and Testing of Pitch-Adapting (Flexible) Composite Propellers[END_REF][START_REF] Liu | Utilization of Bend-Twist Coupling for Performance Enhancement of Composite Marine Propellers[END_REF], and it was also shown that the use of composites can lead to better hydrodynamic performances compared to conventional marine propellers due the bend-twist coupling property. Indeed, this bend-twist coupling enables the propeller to passively adapt to the incoming flow, and thus to increase performances for a wider range of ship speed (Plusincki et. al, 2007). (He et. al, 2012) were then successful in reducing the level of vibrations transmitted to the propeller shaft.

To achieve better performances, it becomes necessary to focus on the design and optimization of such selfadaptable propellers. [START_REF] Lee | Optimized Design of Composite Propeller[END_REF] established the influence of the layers stacking sequence on the hydroelastic performances of composite propellers, and computed the optimal arrangement of the fibers using a genetic algorithm. [START_REF] Liu | Utilization of Deformation Coupling in Self-Twisting Composite Propellers[END_REF]) also related the layers stacking sequence with the capacity of the composite propeller to passively adapt to the flow. In 2009, Liu and Young performed a systematic design method for composite propellers. A year after, [START_REF] Blasques | Hydroelastic Analysis and Optimization of a Composite Marine Propeller[END_REF] realized the optimization of a full-scale prototype and included the structural resistance as a criterion. More recently, (Lee et. al, 2015) optimized the stacking sequence for the investigation of the hydroelastic behavior of the propeller.

With the recent development and validation of effective fluid-structure methods, it is now possible to perform coupled CFD-CSD computations to validate a composite propeller, possibly designed using an optimization process, with consideration for the composite and plies layup and the geometry of the propeller. Structural resistance, cavitation and hydroelastic instabilities can be assessed using common staggered procedures between the codes. This work is the first step of validation of fluid-structure computations which will help the design of composite marine propellers. First, the hydrofoil properties are presented. Then, the governing equations used for the fluid and the structure are given, followed by the numerical models and the set-up used in simulations. The results show the validation against experiments, and then the hydroelastic response is investigated.

EXPERIMENTAL MODEL

The numerical model is taken from the experiments realized by Zarruk et. al in 2014, consisting of a composite and an isotropic hydrofoils tested under the same range of operating conditions in a water tunnel. The hydrofoils are cantilevered at the root and has a free tip that allows the structure to behave like a propeller blade. Details of the dimensions, manufacture, material and structural properties of the hydrofoils, as well as the experimental setup and instrumentation, can be found in (Zarruk et. Al, 2014), [START_REF] Cohen | Numerical Modelling and Experimental Determination of the Dynamic Behaviour of Composite Structures[END_REF] and (Philipps et. Al, 2015). Only the most relevant features of the hydrofoil are reported below.

The isotropic and composite hydrofoils have an identical geometry consisting of a trapezoidal planform with a 300 mm span, and a varying chord in the spanwise direction (Figure 1 and Figure 2). The basechord is 120 mm and the tipchord is 60 mm. The section is a NACA0009 with a modified thicker trailing edge to allow the manufacture of the hydrofoils. To clamp the hydrofoil on the water tunnel wall, the root of the hydrofoil is prolonged with a 20 mm metallic plate, pictured in blue on Figure 2,(A). Two additional metallic plates are used to hold both sides of the hydrofoil in a mounting flange, and a fairing disk is then added to pinch the plates (refer to Figure 1). The incidence of the hydrofoil within the water tunnel can still be varied using this arrangement, and the effective wet span in still 300 mm.

The isotropic hydrofoil was made of aluminum, while the composite one consists of a hybrid design shown in Figure 2 (B). First a thick sandwich glass mat fabric was used at mid-chord to help resin distribution during the RTM (Resin Transfert Molding) process. Then, several plies consisting of unidirectional carbon fiber fabric and biaxial E-glass were stacked on the mat for structural properties. Last, a basket weave E-glass fabric was used as the outer ply to aid surface quality. Since the thickness and the chord are not constant along the span, some plies need to be progressively reduced and dropped relative to the initial stacking sequence shown on Figure 2 (B). This is done similarly to the layers of an onion: the central mat fabric is kept constant along the span, as well as the basket weave E-glass outer ply, and the internal plies the closest to the center are dropped. All fibers are oriented along their longitudinal direction, i.e the 0º-axis. Table 1 sums up the mechanical properties of the different plies for the composite hydrofoil.

[g/cm 3 ] E11 [Mpa] E22 [MPa] υ12 G12 [MPa] G23[MPa] G13[MPa] Orientation [º]

NUMERICAL METHODS

In this section we introduce the governing equations for the fluid-structure problem and the coupling between the two physics, along with the numerical methods used. The numerical set-up is then detailed.

Governing Equations

In this problem we consider a structure Ωs surrounded by a fluid Ωf with an interface Γ, Figure 3. The fluid is characterized by its velocity field v(x,t) and pressure p, and the solid by its displacement field u(x,t), deformation field ε(x,t) and stress field σ(x,t).

Governing Equations for the Fluid

The fluid is modeled using Reynolds-Averaged Navier-Stokes equations (RANS), where the incompressible Navier Stokes equations can be written as (1), with Vi the local velocity and P the pressure.

∂ V i ∂ x j =0 ρ ( ∂ V i ∂ t + ∂ V i ∂ x j V j ) = ρ f i -∂ P ∂ x i +ν ∂ ²V i ∂ x i ∂ x j -ρ ∂ ∂ x j (V i ' V j ' )
(1)

Governing Equations for the Structure

The structure behavior is governed by the fundamental law of the dynamics, the displacement-deformations relation (valid in the case of small deformations), and in the case of composite structures, the anisotropic law of elasticity. These equations are written as (2), with C the stiffness matrix, consisting of 21 independent coefficients.

ρ s ∂ u i 2 ∂t 2 - ∂ σ ij ∂ x j =ρ s g i ϵ ij = 1 2 ( ∂ u i ∂ x j + ∂ u j ∂ x i )
(2) ϵ=C σ

Coupling Equations

Two types of coupling equations (3) are considered: A kinematic condition enforcing the continuity of the velocities in the fluid and in the solid at their interface, and a dynamic condition for the force balance at the interface.

v i = ∂ u i ∂t σ i j ⋅n j =(-p δ i j + μ ( ∂ v i ∂ x j + ∂ v j ∂ x i ))⋅n j (3)

Numerical Set-up and Boundary Conditions

This part describes the fluid and structure models, along with the coupling parameters chosen for the simulations.

The same fluid model can be used both for the isotropic and composite hydrofoils, however, two different structure models have to be set-up.

Fluid Numerical Model

The fluid flow is obtained using the commercial Finite Volume CFD software Starccm+, based on 3D RANS equations. The fluid domain have the same square section and dimensions than the experimental water tunnel (0.6 m high and 2.6 m long). The incoming flow is considered uniform with a reference velocity U0. A symmetry condition (v.n=0) is imposed on the water tunnel walls, whereas a no slip condition (v=0) is set at the hydrofoil's surface (refer to Figure 4). The later can be modified due to the displacements imposed from the structure code.

The numerical clamping is done by locating the base of the hydrofoil on one of the water tunnel wall. The SST kω turbulence model (Menter, 1994) is used in the simulation to accurately predict the boundary layer detachment along the hydrofoil. The domain is discretized using polyhedral elements, with a structured mesh close to the wall to capture the boundary layer. It is particularly refined at the leading and trailing edges in the near wall region and in the wake. An illustration of the mesh is presented on Figure 5 for a section located at the root of the hydrofoil. The boundary layer mesh was performed in order to satisfy a Log Law representation, i.e. 50<y + <70 using a structured grid in the wall vicinity, see Figure 5. Unsteady simulations are considered using a second order backward Euler scheme. The time step is set to satisfy CFL=U0t/xMEAN=1, where xMEAN is the mean cell size in the chordwise direction at the hydrofoil's surface.

Structure Numerical Model

The structure model is setup both for the isotropic and composite hydrofoils. The Finite Element Method is used with an explicit scheme to solve the anisotropic equations and the fundamental relation of the dynamic for the displacements, stresses and deformations fields. Indeed, the time step is identical for the fluid and solid solvers, and is small enough to enable the use of an explicit solver for the structure. This leads to a faster resolution of the structure equations.

To account for the clamping disk, the geometry is extended by 20 mm at the base of the hydrofoil (see Figures 1 and2A), and all the nodes on this additional part are restrained along all directions. The clamping disk is then defined as stainless steel, and solid continuum elements (C3D8) are chosen for the 3D mesh. The numerical model is common between the isotropic and composite hydrofoils, and differs when defining the material properties and the mesh of the hydrofoil. The isotropic hydrofoil is fully modeled with aluminum, and solid continuum C3D8 elements are used for the mesh.

For the composite hydrofoil, these elements cannot account for the stacking sequence of plies described previously and illustrated on Figure 2 (B). Thus the inner and central part of the hydrofoil corresponding to the isotropic sandwich glass mat is modeled classically, using C3D8 solid continuum elements, and a skin with S4R shell elements is added on top and below it (see Figure 6A). The plies are described in the shell elements properties, reported in Table 1. Moreover, as the thickness varies along the span, so does the number of plies: One then needs to cut the hydrofoil into sections, each one of them with a specific number of plies. Figure 6 (A) shows the different sections used to account for the variable number of plies. Fluid and structure meshes on the hydrofoil

Coupling Methodology

In the fluid-structure coupled computation, data is transferred from the fluid to the structure solver to model the exchange of energy between the two physical domains, see Figure 7 (A). In this work, the coupling is realized using the Co-simulation Director Engine (CSE) already implemented in Abaqus to interact with thirdparty codes such as Starccm+, see Figure 7 (B).

Two types of calculation are performed, based on different coupling algorithms. A non coupled procedure is first set. It uses only one converged fluid solution, then exports the pressure field to get the structure solution. Tightly coupled fluid structure computations are also performed. They are based on a two way coupling between the fluid and structure codes, where staggered iterations within one coupling time step help for the energy conservation at the interface, as illustrated on Figure 7 (C). To advance in time, the displacements and pressure must indeed reach a state of convergence between the values (u,p) k t-1 and (u,p) k+1 t-1, or be forced to move forward in time after ten exchanges, where k corresponds to the sub-iteration number. Moreover, the fluid resolution scheme is implicit, so for each data exchange, the pressure p k t is determined using several fluid internal iterations to reach convergence.

To increase the stability of the coupled calculation, the displacement solution imported in the fluid solver is under-relaxed to progressively implement the current value during the coupling time. The first two data exchanges are not under-relaxed, then the relaxation factor is optimally adjusted between 0.2 and 0.5. Moreover, the coupling between the two solver is also progressively implemented for stability issues. Figure 8 shows the progressive implementation starting with a fluid steady calculation for the case of an infinitely rigid hydrofoil (no coupling -in brown on Figure 8) until convergence. An unsteady calculation is then performed on the same case until a pseudo-convergence is reached (in blue). The results act as an initial condition for the coupled calculation, increasing the initial stability.

RESULTS

This part details the results of the different simulations, including the validation of the fluid and structure models separately, and the results from non coupled and coupled simulations.

Fluid Model Validation

The fluid model is first checked for the case of infinitely rigid hydrofoil (no coupling).

Mesh Convergence Study

A mesh convergence study on the fluid domain is performed for 8° angle of attack. Seven meshes are considered, with numbers of elements ranging from 217,810 to 4,004,628, where the mesh topology and boundary layer mesh is kept, see section 4.1 and Table 2. The convergence of hydrodynamic coefficients are presented on Figure 9. It shows that the solution is converged for the lift coefficient starting from 543,535 elements.

Structure Model Validation

To validate the structure model, three criterias are checked. First, the weight of the hydrofoil is numerically adjusted to match the experimental weight by tuning the density coefficient of the materials. Then, natural bending and twisting frequencies are adjusted using Young's Modulus E22 and the Shear Modulus G12 of the structural materials (unidirectional carbon and E-glass weave). For the isotropic hydrofoil, only experimental natural bending frequency is provided, and numerically matched. For the composite hydrofoil, the validation of the above three criteria is presented in Table 3. he maximum tip displacement of the hydrofoil under hydrodynamic loading is then determined for the non coupled and the coupled calculations, and compared against experimental measurements. Results are reported in Table 4 for the isotropic hydrofoil, and in Table 5 for the composite hydrofoil. There is a general good agreement between the present computations and measurements. For the aluminum hydrofoil, the computed drag coefficient is above the experimental value. There is also an under prediction of the tip displacement. The results for the coupled calculation, which will be presented in the final version of the paper, are expected to get closer from the experiments due to the effect of fluid structure interaction, i.e. the coupling between twist deformation and lift coefficient. For the composite hydrofoil, the tip displacement is also under predicted both at the leading and trailing edges. However, the twist predicted on the hydrofoil presents good agreement with the measurements. As the twist is greater for the composite hydrofoil, due to its increased flexibility, fluid structure interaction effects are expected to be more important in the coupled calculation than for the aluminum. As a consequence, it is expected that the coupled results will show a significant increase of lift force together with an increase of bending and twist deformation, which will increases the tip displacement as well. Comment to reviewer: Strong coupling computations are currently being processed, and results will be reported in Table 4, replacing the X, in the paper final version. In this work, the flow over flexible NACA009 3D hydrofoils is simulated using a finite volume/finite element coupling between the CFD solver Starccm+ and the CSD solver Abaqus. Two hydrofoils of the same dimensions and trapezoidal planform geometry, but of different materials, were simulated in a cantilevered arrangement to assess their hydroelastic behavior, and results were compared with experiments from Zarruk et al. 2014. The first hydrofoil was made of aluminum, and the second of carbon-fiber reinforced plastic (CFRP) with fibers oriented at 0º. Fluid models and structure models were setup, preliminary validated, and then computed using non coupled and coupled procedures. The results were validated on the lift and drag coefficients, and the maximum tip displacement of the hydrofoil. It shows overall very good agreement with the reference study. Moreover, the numerical study confirmed that the composite hydrofoil has a higher twist deformation than the aluminum hydrofoil, due to the 0° fiber orientation along the span, whereas it seems more rigid if no coupling is considered. Natural bending and twisting frequencies are adjusted using Young's Modulus E22 and the Shear Modulus G12 of the structural materials, and this might explain the increased stiffness of the composite hydrofoil. However, the larger twist deformation will tend to amplify the coupling effect induced by the local increase of the angle of attack along the span, which interacts with the lift force. The coupled results for steady and unsteady flows will be presented in the final version of the paper.

Figure 1 .

 1 Figure 1. Schematic illustration of the general mounting arrangement of the hydrofoil for water tunnel tests, from (Cohen et. al, 2014)

Figure 3 .

 3 Figure 3. Theoretical formulation of the fluid-structure interaction problem.

Figure 4 .

 4 Figure 4. Schematic illustration of the fluid domain and boundary conditions.

Figure 5 .

 5 Figure 5. Illustrations of the fluid domain mesh

Figure 6

 6 Figure 6. (A) Model of the composite hydrofoil showing the different sections used to account for the variable number of plies within the hydrofoil, and the mesh elements types, (B)

Figure 7 .Figure 8 .

 78 Figure 7. (A) Principle of a fluid-structure coupling, (B) Numerical Codes and Methods used in this work, (C) Details of the Staggered Implicit Coupling Algorithm

Figure

  Figure 9. Mesh convergence study using seven meshes for the configuration with 8°angle of attack and Re = 600,000

  results obtained for the non coupled and coupled calculations for an angle of attack of 8°. Results of the non coupled calculation are first used to validate the accuracy of the pressure transfer between the fluid and solid solvers, using a Least Squares interpolation method. Figure11shows the pressure computed with the CFD code and the mapped pressure read in the CSD code, relative to the fluid and structure meshes shown on Figure6 (B). It is shown that the pressure gradient at the hydrofoil's surface is correctly interpolated in CSD code, whereas the mesh topology is quite different.

Figure 11 .

 11 Figure 11. Validation of the transfer of pressure between the fluid solver and the structure solver, =8°T

  Figure 12. (A) Displacements along the Y-direction for the aluminum hydrofoil and (B) for the composite hydrofoil, =8°T able 4. Results of the fluid-structure computation on the isotropic (aluminum) hydrofoil, =8° -The values for strong coupling calculations, marked with an X, will be provided in the final version of the paper Aluminum hydrofoil CL CD Max. tip Displacement [mm]

Table 1 . Properties of the Composite Plies

 1

	Carbon E-glass	Basket	Mat
	ρ

Table 3 . Summary of structure model validation

 3

Acknowledgements

The authors thanks DCNS for its financial support through a collaborative grant DCNS-LHEEA. The calculations were performed on a Nehalem cluster at CRIANN, France.


M Atlan

M Gross

B C Forget

T Vitalis

A Rancillac

A K Dunn

Frequency-domain, wide-field laser Doppler in vivo imaging

We present a new instrument, based on a low frame rate (8 Hz) CCD camera used in a heterodyne optical-mixing configuration, that can create wide-field laser Doppler maps. As an illustration, we show results obtained in a mouse brain, in vivo, showing the Doppler signature of blood flow. The instrument is based on the frequency-shifting digital holography scheme.

Imaging blood flow is an essential tool to assess many physiologic processes and diseases. Several instruments have been proposed to produce blood flow maps. These instruments can be grouped in two categories : those which are based on time domain measurements (laser Doppler techniques) and those based on spatial domain measurements (speckle contrast analysis techniques) [1,2]. The first Laser Doppler experiment is attributed to Yeh and Cummins [3], who measured the frequency shift of light scattered by moving particles with a spectrometer.

Recently, Serov et al. proposed [4] and designed [5,6] a full-field Laser Doppler imaging technique based on the acquisition of a sequence of images by a fast CMOS detector. This method requires on the one hand a frame rate high enough to set the cut-off frequency of the restituable spectrum higher than the ones induced by typical blood flow, eg. [6] ∼ 20 kHz and, on the other hand, a sequence of images large enough to reach a 50 Hz spectral resolution [5]. Laser speckle contrast analysis [7] has been proposed as an alternative to scanning laser Doppler imaging as a full-field imaging technique. In this technique, the first order spatial statistics of the speckle pattern are used to form blood flow images. It has been successfully used in applications including noninvasive skin capillary blood flow imaging [8] and functional brain imaging [1]. Laser Doppler and speckle contrast analysis turn out to be based on the same underlying phenomenon [2,[9][10][11]. A new experimental approach is proposed here, consisting of a frequency domain measurement onto a parallel detector. This measurement is based on an optical mixing scheme.

The setup consists in a heterodyne measurement of scattered (object) optical field E O . The object is illuminated over a ∼ 1 × 1 cm 2 area (the region to be imaged). The backscattered laser wave is mixed with a frequency-shifted local oscillator (LO) E LO field. Laser light is provided by a 658 nm (angular frequency ω L ), 80 mW CW Mitsubishi ML120G21 single-mode diode. Two Bragg cells (Acousto-optic Modulators, AOMs), are used as beam frequency shifters (see fig. 1). The LO angular frequency ω LO is shifted by the difference ∆ω AOM in driving frequencies of the set of Bragg cells : ω LO = ω L + ∆ω AOM . A 1.3 Megapixel PCO PixelFly camera (framerate ω S , square pixels, pixel size 6.7 µm) is used to record a sequence of m images (namely I k , where k = 1, . . . , m) of the interference pattern resulting from optical mixing between the LO and object fields. The detection plane is set 50 cm away from the object plane. A circular pupil set in front of the object limits the field of view. Since the spatial field of interest in the object is ≈ 1 cm 2 , it matches the natural heterodyne parallel detection image field, avoiding the need for any lens in the object path. The optical configuration is a typical holographic setup, called lensless Fourier (off-axis) holography setup [12,13], in which the point source of the spherical reference wave (LO focal point) is located in the plane of the object. The lensless Fourier setup is used to avoid the need of a general Fresnel reconstruction algorithm [14] to assess the spatial distribution of the scattered field in the object plane. In this case, the reconstruction algorithm consists of only one spatial fast Fourier transform [13] (FFT). The digital hologram is recorded according to the frequency-shifting method introduced in the heterodyne holography technique, used in several imaging and detection schemes [15][16][17][18] and which turns out to be sensitive enough to discriminate Doppler-shifted photons, according to their frequency, from background light [16,18]. The frequency-shifting method was chosen as an accurate [15] phase-shifting [19] detection method. To create a map of the ∆ω angular frequency component (corresponding to the ∆f = ∆ω/(2π) frequency, in Hz) of the scattered field in the image plane, the LO frequency is set to ω LO = ω L + ∆ω + ω S /n, where the integer n satisfies n ≥ 2. The power S(ω = ω L + ∆ω) of the ∆ω frequency component of the measured object field is computed in parallel on the pixels of the image, according to the n-phase demodulation method [16,17]. It consists in calculating the squared amplitude of the ω S /n AC component of intensity fluctuations I k by a time-domain FFT:

Laser AOM1 AO M2 camera cube M 1 M 2 M 3 M 4
S(ω) = A | n k=1 I k exp (-2j(k -1)π/n)|
2 , where A is a constant and j 2 = -1.

In the reported experiment, the cranium of the mouse is illuminated with an incidence angle α ≃ 60 • . The angle between the object and LO arms is set to θ ≃ 1 • . This experiment was designed to determine whether the instrument would allow detection of the cerebrocortical blood flow.

Animal procedures were conducted in strict compliance with approved institutional protocols and in accordance with the provisions for animal care and use described in the European Communities Council directive of 24 th November 1986 (86-16-09/EEC). Five C57/Bl6 mice (Charles River) aged P25-P30 were anaesthetised with a mixture of Xylazine-Ketamine (1 mg/kg IP, 10 mg/kg IP). Animals were positioned on a stereotaxic frame (World precision Instruments) to ensure stability of the preparation. Cranial skin and subcutaneous tissue were excised linearly over the sagittal suture and cortical bones were preserved. To assess bioflow [20], the root mean square (RMS) frequency shift ∆ω 2 1/2 is calculated with respect to the S(ω) distribution measured at 80 linearly spaced frequency points in the 0 → 2528 Hz range (spacing interval : 32 Hz). The camera framerate was set to 8 Hz. A sequence of m = 32 images sampled during 32/8 = 4 seconds was taken to build each S(ω = ω L +∆ω) map at one given ∆ω frequency point, setting ∆ω AOM = ∆ω+ω S /n, where n = 4. To increase SNR, S(ω) was averaged by means of an accumulation of sliding n-phase demodulation along the m image sequence. The whole measurement took 5 minutes and 20 seconds (= 32 images × 80 frequency shifts × 1/8 second exposure time). Two S(ω) spatial maps at a low (∆f = 64 Hz) and higher (∆f = 2240 Hz) frequency shift and a spatial map of the RMS the frequency shift (in Hz) are represented on fig. 2. Three spectra, averaged over 5 × 5 pixels re- gions of interest (ROI, outlined in fig. 2) are displayed in fig. 3. The larger blood vessels (eg. ROI 3) induce a broader spectrum than the smaller ones (ROI 2). In the parenchyma (ROI 1), spectra are narrower, suggesting a reduced blood perfusion. Hence the contrast reversal between the two spectral maps in fig. 2. An additional in vitro calibration experiment was done to assess the order of magnitude of the measured RMS the frequency shift for various flow rates. The result is shown in fig. 4; it demonstrates an approximately linear relationship between the average speed (from 0 to 750 µm.s -1 ) and the measured frequency shift. In conclusion, this scheme based on heterodyne optical mixing onto a parallel detector is particulary suited to wide-field laser Doppler imaging. An arbitrary frequency component of the object field can be acquired at once. The available range of frequency shifts at which measurements can be made is far more extended than in any of the concurrent techniques, since the measurement is directly performed in the frequency domain. The actual limit lies in the system used to shift the LO frequency. As an indication, the RF-driven AOMs combination can shift this frequency up to 10 Mhz. The measurement is sensitive : thanks to the off-axis holographic setup, the low frequency noise components (homodyne intensity terms, resulting from laser instabilities and scattered light self interference) within the camera bandwidth are efficiently filtered-out [16,17]. The Doppler map is not restituted over the total number of pixels of the camera. Only up to 1/4 of the pixels of the frame can be used to describe the image, as a result of the off-axis holographic scheme. Nevertheless, both temporal and spatial resolution are potentially high. This combination is enabled by a measurement performed at one frequency point at a time. A faster camera could be used to improve temporal resolution and still keep a spectral resolution acceptable for this measurement, since the spectral resolution is defined by the heterodyne bandwidth (inverse of the acquisition time). A detection array with a large number of pixels should be used to perform an efficient wide-field and high spatial resolution detection. Due to the data transfer rate limit encountered, the balance between the number of pixels and the frame rate of the detector should be guided by the application needs in terms of temporal, spectral and spatial resolution.

Contact : atlan@lkb.ens.fr

FIG. 1 :

 1 FIG. 1: Optical configuration (general view, upper fig. Inset, lower fig.). AOMs : Bragg cells. BS : beam splitter. BE : beam expander. M1 to M6 : mirors. EO : Object (scattered) field. ELO : reference (LO) field.

FIG. 2 :

 2 FIG. 2: Power maps of the ∆f = 64 Hz and ∆f = 2240 Hz frequency components of the scattered field (top, log. scale) and perfusion map, i.e. RMS frequency shift, in Hertz (bottom). The image shows a dorsal view of the mouse cranium (anterior on the left, posterior on the right). The superficial dorsal venous system and some of the superficial cerebral arteries are visible.

  FIG. 3: Spectra calculated in the three 5 × 5 pixels ROI outlined in fig. 2 (log. scale).

FIG. 4 :

 4 FIG. 4: In vitro measurement of the light scattered by a 5 10 -3 vol. frac. suspension of 1 µm latex beads flowing through a transparent tube of 580 µm of diameter. Average RMS frequency shift vs. average speed of flow.


Adrien Van Den Bossche

Thierry Val
email: val@iut-blagnac.fr

Eric Campo
email: campo@iut-blagnac.fr

Proposition of a full deterministic medium access method for wireless network in a robotic application

Keywords: component; QoS, Wireless, LP-WPAN, ZigBee, IEEE 802.15.4, MAC, Sensor networks, Mesh networks, Mobile adhoc networks I

Today, many network applications require shorter react time. Robotic field is an excellent example of these needs: robot react time has a direct effect on its task's complexity. Here, we propose a full deterministic medium access method for a wireless robotic application. This contribution is based on some low-power wireless personal area networks, like ZigBee standard. Indeed, ZigBee has identified limits with Quality of Service due to non-determinist medium access and probable collisions during medium reservation requests. In this paper, two major improvements are proposed: an efficient polling of the star nodes and a temporal deterministic distribution of peer-to-peer messages. This new MAC protocol with no collision offers some QoS faculties.

INTRODUCTION

Today, networks and telecoms are widely used. In production, space technologies, or even at home, networks and moreover wireless networks are considered to become essential. Wireless technologies eliminate expensive, heavy, not aesthetic cables, which are not easy to install or to use. Using a wireless connection is sometimes a luxury, but it may be necessary in many cases of moving devices, like car tire sensors, etc. All these points encourage research and industrial to develop technologies and products in this domain.

Nowadays in telecoms and networks fields, developers are essentially trying to make the products cheaper and easier to use. Baudrate is constantly increasing and there is currently a high level of performance and compatibility between technologies which are globally mature. Unfortunately, there is a technical problem which is not yet solved: temporal performance on the delivery of the network messages. In the 80's, Internet network applications did not need good temporal performance because there were no real-time applications over the network. Today's applications require higher temporal performances: audio/video streaming for home applications, sensors networks [START_REF] Callaway | Wireless Sensor Networks: Architectures and Protocols[END_REF], medical assistance [START_REF] Golmie | Performance evaluation of low rate WPANs for medical applications[END_REF], etc. are widely used and imply temporal certitudes; medium access needs to be determinist and a maximum latency has to be kept within identified limits. In many networks, errors are generally corrected by protocols with an acknowledge/retransmit mechanism; for many real-time applications, this method can't be used because of the latency introduced by retransmissions [START_REF] Van Den Bossche | Metrologie pour l'analyse comparative des performances temporelles des liens Bluetooth[END_REF].

Another point that can be noticed in this introduction is the interest to regroup network and embedded sciences. Today's wireless devices are portable; because of it, energy saving becomes a major goal in the development of an electronic device. This point is crucial and developers have to ask themselves about conventional methods for solving classical problems. For example unnecessary retransmissions, padding, frame collisions or permanent listening of an unused medium are energy losses which must be limited. IEEE recently worked on Low Power Wireless Personal Area Network (LP-WPAN) and some new technologies like IEEE 802.15.4/ZigBee can solve many problems on simple communicating embedded devices. This work plans to improve network protocols for optimal energy saving.

II. OUR WIRELESS ROBOTIC APPLICATION

A. General presentation and objectives

Our application takes place in mobile robotic field. Mobile robots are obviously embedded devices. They can communicate with a central entity of command which can be mobile/fixed or with other robots. In that case, robots are cooperating. It has two goals: first, eliminating many wires in mobile robots which are equipped with sensors. Research prototypes robots will be easier to repair and to upgrade. Moreover, robots will be lighter which can be fundamental for unearthed vehicles. The other goal of the application is the creation of communication channels between robots by this internal network. By this way, cooperating robots can listen to messages of other robots and use this listening for learning about their environment. This capacity is interesting for cooperating robots: regulate pursuit speed, mapping environment, etc. For example, if one of the robots breaks down because of an external physical problem (due to a fall, excessive heat, clash, etc.), others robots can intercept the message of the sensor and react to help or protect themselves.

B. Network proposed

As an illustration of this work, our application plans to create two types of communication channels as illustrated in fig. 1: between a group of sensors/actuators and the robot command unit, and between robot command units.

For the first type of communication, the network will naturally be configured as a star topology with a unique star in each mobile robot. Each star is composed of a coordinator located in the command unit of the robot and of network nodes for each sensor. Each node establishes communications with the coordinator via star links and cannot communicate directly with another node without going through the coordinator.

Mobile robot, i.e. command device, can also communicate with other mobile devices or fixed network infrastructure by establishing peer-to-peer network links. Only coordinators can create peer-to-peer links, nodes cannot. The application may also support coordinators without nodes, if slave nodes were all turned off or if there was no node present for this coordinator in the area. To test this network, we recently designed several network adapters based on the Freescale semiconductor MC13192 chip [START_REF] Semiconductors | MC13192 2.4 GHz Low Power Transceiver for the IEEE 802.15.4 Standard: Reference Manual[END_REF]. These devices are placed on some ER1 robots, by Evolution Robotics [START_REF]Official website and ER1 overview[END_REF].

III. OVERVIEW OF IEEE 802.15.4 STANDARD

In the last years, many IEEE's works study Low Power Wireless Personal Area Networks (LP-WPAN) such as IEEE 802. 15.4 [6]. This standard proposes an original two-layer protocol stack (physical-layer and data link-layer) for low power transceivers and low baudrate communications between embedded devices. Energy saving has been optimized by using innovative concepts.

A. Nodes types

An IEEE 802.15.4 network is built with several types of devices: simple nodes, routers, star coordinators or network coordinators. 802.15.4 standard describes two versions of the protocol stack. On the one hand there is a Full Function Device (FFD) stack which includes all functionalities proposed by the standard. A FFD can achieve all functionalities of the network: terminal node, routing and coordinating. On the other hand there is a lighter version of the stack named Reduced Function Device (RFD). RFDs can be only a terminal node of the network. There are few levels of reduction of the stack: with or without beacon mode ( §III.C.2), with or without crypt functions, etc. Generally, FFDs are main powered and RFDs are embedded.

B. Physical layer

IEEE 802.15.4 physical layer has been designed for maximum energy saving:

• DSSS (Direct Sequence Spread Spectrum) O-QPSK (Orthogonal-Quadrature Phase Shift Keying) modulation: excellent noise immunity so transceiver uses low power over the air,

• good receiver selectivity (-95dBm),

• low baudrate: a high noise immunity, few transmission errors so few treatments for higher levels (link and network); using simple processors (8bit) is possible.

All 802.15.4 nodes use the same radio channel. Protocols are optimized for short and periodical data transfers: nodes are sleeping most of the time. In this IEEE standard, this mode is called doze mode. It allows ultra low power consumption (40µA) and nodes can become operational in a very short time (330µs) [START_REF] Zarader | An overview of ZigBee/802.15.4 Freescale Solution[END_REF]. In doze mode, all radio functionalities are powered down, so it is not possible to receive messages over the air. Dozing devices must establish together a future waking moment before switching doze mode (synchronous wake up), but sleeping devices can also wake up on a local event (asynchronous wake up), for example on a sensor detection. In a conventional use of the network, devices sleep more than 90% of the time. Scheme is: devices are sleeping, then wake up, transmit/receive data, then at last go to sleep again (see §III.C.1). Only coordinators stay awake for buffering network messages of sleeping terminals.

C. Data link layer

Data link layer has two objectives: medium access control and detecting/correcting transmission errors. 802.15.4 also proposes two link level topologies.

1) Topologies: 802.15.4 proposes two topologies: peerto-peer and star.

• Peer-to-peer topology makes possible direct data transfers between FFD. Establishing peer-to-peer links is only possible if devices are in common radio range and use the same radio channel. Medium access is done by using CSMA/CA protocol without RTS/CTS mechanism ( §III.C.2). Peer-to-peer topology makes ad-hoc networks [START_REF] Murty | Ad Hoc Wireless Networks: Architectures and protocols[END_REF] possible.

• Star topology is possible when a coordinator is present on the channel. The coordinator is a FFD which manages medium access of "slave" nodes. Coordinator's "slave" node can be either a FFD or a RFD. All data transfers are going through the coordinator, making doze mode possible -messages are buffered by the coordinator. This network topology allows high energy saving thanks to an optimal sleeping period distribution between embedded devices. Star coordinator sends beacons frames on demand or periodically; periodically beacons frames make beacon tracking possible, a special mode which consists of a RFD in a slept waiting for the next beacon frame, then wakes up, asks the star coordinator for pending data, transmits/receives, and goes to sleep again.

2) Medium access methods: IEEE 802.15.4 proposes two medium access methods:

• A method based on CSMA/CA contention method,

• A contention free method which allows bandwidth reservation.

The contention method can be used in either peer-to-peer or star topologies but contention free method can be used only in star networks. To use contention free method, star coordinator has to attribute Guaranteed Time Slots (GTS) to its nodes, others nodes cannot transmit during these allocated periods. Thanks to contention free medium access, medium reservation is possible and some Quality of Service (QoS) properties like bandwidth reservation or latency guaranties can be applied. This contention free access is only possible in a star beaconed network. Beacon frame contains information which indicates the structure of inter-beacons period named superframe. IEEE 802.15.4 superframe has the following structure:

• At first, star coordinator sends a beacon frame which indicates durations of superframe and sleeping interval, data pending, network and coordinator addresses, cryptographic properties, size of the Contention Access Period (CAP) and size of the optional Contention Free Period (CFP),

• Coordinator and nodes send/receive frames in the CAP using CSMA/CA protocol (GTS requests are done in the CAP),

• Coordinator or nodes use GTS in the optional CFP,

• Before next beacon frame, the optional sleeping period can be imposed by the star coordinator for energy saving. All network nodes go in doze mode (see §III.B) and local timer will wake up each device just before next beacon.

D. ZigBee: higher layers

ZigBee is a standard of the ZigBee Alliance [START_REF]ZigBee Alliance[END_REF] which uses IEEE 802.15.4 standard for physical and data-link layers. ZigBee proposes an original protocol suite for 3 to 7 protocol stack levels [START_REF]ZigBee Specifications[END_REF]. With its network layer, ZigBee extends network range and can obtain wide zones, using mesh architecture for example with cluster-tree topology [START_REF][END_REF]. Thousands of devices can be linked by this way. With its profile system, ZigBee provides a high-level of compatibility between manufacturers. Profiles concept was created for IrDA, and has been widely approved with Bluetooth systems [START_REF] Palowireless | Bluetooth profiles[END_REF].

E. IEEE 802.15.4 in our application

According to what was mentioned above, IEEE 802.15.4 is an interesting technology for our application. Topologies adapted to our problem, low power, data encryption, sufficient baudrate and, of course, restricted mobility is possible.

1) Topologies: Topology possibilities of 802.15.4 are adapted to our needs: a robot will be equipped with a star network, robot command and control device will be the star coordinator (FFD) and sensors/actuators will be simple network nodes (RFD). Robots can also communicate with other robots via peer-to-peer links using ZigBee mesh topology.

2) Mobility: In this application, mobility between different robots is required but will be limited inside the robot: sensors/actuators will not move within the robot or only with a limited mobility, for example a robot arm sensor. Cooperating robots will be coordinated, so mobility is also limited for peerto-peer network links. Thus network has certain topology stability. 802.15.4 supports "limited mobility" (Access and association delays, time to live and route validity periods)

3) Energy saving: Embedded characteristic of the robot makes this point crucial. Other IEEE network technologies like 802.11 (WiFi), 802.15.1 (Bluetooth) are not advanced enough in this respect.

4) Baudrate and latency:

The application does not need high baudrate but messages must pass through the network with guarantee on delivery time (latency). Network carries crucial sensors and actuators data (obstacle or shock sensor, "stop" motor order, etc.) which can not tolerate delays and must be forwarded in time. Several delay guarantee levels must be defined over the network but medium access without delay guaranty also has to be preserved because all nodes do not need time guarantees. Moreover, the application generates a little traffic (few kilo-octets per second) but it has a very consistent shape. Sensors raw messages are small (generally few octets) and temporal needs are so strict that data cannot be buffered in order to make full network frames. Because of that point, efficient data are very small and frame headers must also be as small as possible. Efficient baudrate will be much reduced as well. The network technology must be the simplest on that point and protocol must be light. 802.15.4 has been conceived in order to present simplest protocol that can be implemented on small 8-bit microcontrollers, that is again an advantage for the application.

F. Limits of IEEE 802.15.4 in our application

802.15.4 is an interesting wireless technology but, unfortunately, it is not perfect. For example, in a star network present standard cannot certify deterministic medium access. As described in §III.C.2, it is possible to reserve GTS but GTS requests frames are sent in the CAP, where frames may generate collisions. This gap is embarrassing, it is impossible to guaranty transmission delays. Even worse: CSMA/CA protocol used in the CAP is based on random message delaying to avoid collisions, which is not adapted for real-time applications. 802.15.4 has another gap: GTS cannot be used in peer-topeer links. Again, it is not possible to reserve bandwidth, guarantee transmission delay is not possible. These gaps put a stop to QoS exchanges over the network.

The next section presents our proposition on 802.15.4 to eliminate these two disadvantages thanks to a medium access method modification.

IV. NEW MEDIUM ACCESS MEDIUM PROPOSED

Works presented in this section aim to solve 802.15.4 gaps presented in last section: possibility to make determinist GTS requests i.e. without collisions and to avoid collisions in peerto-peer network links.

A. GTS requests without collisions

802.15.4 standard makes possible nodes to request GTS to the star coordinator in the CAP, so these requests may cause collisions. To solve this problem, we propose here to reverse the direction of the request. For some important nodes, star coordinator can ask the associated node for possible GTS allocation, rather than keep node doing its GTS request; requests can be done by the coordinator in a polling. Coordinator must have pre acquired knowledge of its nodes and their needs. In our application this assumption is realistic. This new method is an improvement which permits to guarantee a maximum delay for GTS allocation. Withoutpriority nodes can still make GTS requests in the CAP and vital nodes can also accelerate their requests by making requests in the CAP.

1) The Network association requests and GTS requests:

In the life of an IEEE 802.15.4 network, there are two critical temporal phases: network association and GTS requests. In present standard, network association i.e. coordinator association is achieved by a passive listening of the radio medium, searching for beacons frames. All channels must be scanned to find a suitable coordinator; each channel listening period must be calculated according to the beacon transmitting frequency, which must be declared in the node configuration [START_REF] Neugeabauer | A new Beacon Order Adaptation Algorithm for IEEE[END_REF]. Once network association is completed, a network node can transmit/receive data, goes in doze mode or request for GTS. These two phases of an IEEE 802.15.4 network life (network association and GTS Requests) are very close to each others if we are searching for temporal guarantees: it may be interesting to obtain a definite response to these requests. Moreover, a definite response is the best to consider energy saving.

2) Pre acquired knowledge of the neighborhood nodes:

Star coordinator must have pre acquired knowledge of the neighborhood nodes and also the needs of the sensors for making a perfect polling. For our application, this assumption is not restricting because it is the command unit of the robot which organizes the network. This central entity must know all the sensors it has to take care of them. Pre acquired knowledge can be a list of MAC addresses, with a precise description of this node's needs: minimal and maximum baudrate, maximum latency, etc. This knowledge base can be downloaded in the coordinator's memory or created by automatic learning. Controlling first network association can be done by pushing manually synchronization buttons like in DECT systems.

3) Polling optimization: Once pre acquired knowledge primitives are implemented, polling has to be optimized: simple linear polling, scheduling with priorities, etc. Several algorithms were developed to solve these classic problems. Each application has its adequate algorithm. Each network star can have its own scheduling.

B. Peer-to-peer links without collisions

In the previous section, we saw how to solve the deterministic medium access gap of IEEE 802.15.4 standard for GTS requests in a star network topology. In this section, we propose a solution to avoid collisions in a multiple peer-to-peer links area. 802.15.4 standard does not propose a GTS mechanism for peer-to-peer links. In today standard, FFDs can communicate outside the star network with peer-to-peer links using CSMA/CA protocol but because of CSMA/CA employment, present standard guarantees neither baudrate nor transmission delays.

1) GTS between star coordinators: Our solution uses the concept of star's GTS. In a star network, each node is associated to a coordinator; this coordinator decides of the medium temporal distribution between nodes. In peer-to-peer links, there is no coordinator; All nodes which are capable of communicating directly are FFD. In this case, it is possible for two FFD to communicate by using the star mode rather than peer-to-peer mode. If a coordinator A wants to send a message to a coordinator B, A must be previously associated to B and requests for a GTS of B's star. In the next superframe, A can send its message to B without collision. GTS temporal position is sent in B's beacon so A must listen to that beacon. Consequently, it is necessary for each coordinator to listen to the next beacon after its GTS request. If a third coordinator C is in the neighborhood, it may not use the medium while B is sending its beacon. Consequently, it is necessary for all network nodes to keep medium unused while beacons messages are being sent. A hard synchronization is required for our temporal organization of the network.

2) An equal temporal repartition of the beacons: As a starting point, our solution needs each coordinator to listen to every neighbor's beacons. At first, the hypothesis that all coordinators can listen to all others coordinators is mandatory; in that case, there must be no hidden coordinator. For a robotic application, this hypothesis is not so restrictive. However this simplification must be temporary because our contribution could be used again for another application, for example in a wider area. The hidden coordinator problem can be solved by adding a supercoordinator to our network (for a centralized network solution) which identifies hidden beacons or by creating synchronization and beacon relaying messages between coordinators (for a decentralized network solution). The IEEE 802.15.5 task group proposes to avoid beacon collisions of neighbor WPANs to keep independent traffics. It is done by dynamic beacon alignment [START_REF] Dacosta | Dynamic Beacon Alignment in Simultaneously Operating Piconets using the Heart Beat approach[END_REF] which consists in hidden coordinator's beacons retransmission.

By using a wide synchronization over the network, simple nodes must also be synchronized: nodes mustn't send frames during beacon transmitting periods. It is necessary to impose reserved timeslots for beacons and reserved timeslots for nodes. The question is to know how beacons reserved timeslots must be positioned to ensure the best repartition of beacons in the superframe and avoid GTS freezing. This problem may appear if there are too many coordinators in the network area: medium will be saturated by beacons; then the first coordinator can allocate all timeslots for its nodes and medium will be immediately full. A first solution consists in an equal time distribution of beacons, as shown in figure 2.

3) Traffic optimization -coordinator's GTS requests in beacons:

As shown in figure 2, medium temporal distribution is quite condensed. Useless messages have to be limited to optimize bandwidth. In 802.15.4 standard, beacon message size varies with the number of nodes data pending and number of GTS. In a low loaded star, beacon messages are empty. Because beacons are listened to by all coordinators, it may be interesting to use these messages to transport useful information for all coordinators: neighbor coordinators list, received power report, errors report, etc. GTS requests between coordinators can be done via a two beacons request/reply rather than the present standard method, to limit bandwidth loss.

4) Optimization -variable transmission powers:

Figure 2 shows that various star internal administration messages (GTS requests, acknowledgment, etc.) use baudrate and are not useful for other stars. Another optimization for a coordinator consists in allocating GTS to other coordinators GTS after a preliminary study of received power signal of the beacon or the node message: if a first coordinator allocates some GTS to a node that its messages are always received with a low power by a second coordinator, this coordinator can also allocate a GTS in this timeslot, asking its node to transmit in this time slot with a minimal transmitting power. Two coordinators may allocate the same time-slot for two faraway RFD without making collisions. Medium occupation could be denser as shown in figure 3. The work deals with an improvement of actual 802.15.4 standard by adding some QoS features. The proposition consists in a full deterministic medium access method, keeping compatibility with official version of the standard. Two main features are realized: an efficient polling of the nodes by the star coordinator and a temporal deterministic distribution of peer-to-peer messages in a mesh topology. We are validating this protocol to prove its efficiency by adding our contribution in existing 802.15.4 NS2 model [15]. In the same time, we are working on a prototype based on a Freescale development kit to implement our modifications on a couple of real devices. A concrete testing application is also being developed in a robotic prototype [START_REF] L²i Group | Wireless for roBOT (WiBot) project[END_REF] in coordination with the Protocols & Networks Research Team of LIMOS Laboratory, University of Clermont Ferrand, France. Future works are numerous: full validation of the network, prototype testing, improvement of the protocol, solving the hidden coordinator problem.

Figure 1 .

 1 Figure 1. Network topology inside and between robots

Figure 2 .

 2 Figure 2. Temporal repartition for 3 stars

Figure 3 .

 3 Figure 3. Optimized medium use

ACKNOWLEDGMENT

These research works are deal with the industrial collaboration of Cyril Zarader, Freescale Semiconductor, Toulouse, France and LIMOS Laboratory, University of Clermont Ferrand, France.


J Abascal

N Ducros

V Pronina

S Bussod

A Hauptmann

S Arridge

P Douek

F Peyrin

MATERIAL DECOMPOSITION PROBLEM IN SPECTRAL CT: A TRANSFER DEEP LEARNING APPROACH

Keywords: Spectral CT, inverse problem, deep learning, transfer learning

Current model-based variational methods used for solving the nonlinear material decomposition problem in spectral computed tomography rely on prior knowledge of the scanner energy response, but this is generally unknown or spatially varying. We propose a twostep deep transfer learning approach that can learn the energy response of the scanner and its variation across the detector pixels. First, we pretrain U-Net on a large data set assuming ideal data, and, second, we fine-tune the pretrained model using few data corresponding to a non-ideal scenario. We assess it on numerical thorax phantoms that comprise soft tissue, bone and kidneys marked with gadolinium, which are built from the kits19 dataset. We find that the proposed method solves the material decomposition problem without prior knowledge of the scanner energy response. We compare our approach to a regularized Gauss-Newton method and obtain a superior image quality.

INTRODUCTION

The new generation of spectral computed tomography (SCT) scanners include photon-counting detectors, which count single photons and resolve their energy [START_REF] Roessl | K-edge imaging in x-ray computed tomography using multi-bin photon counting detectors[END_REF]. With this extra dimension, SCT allows for material decomposition, which opens up new diagnostic possibilities [START_REF] D P Cormode | Atherosclerotic Plaque Composition: Analysis with Multicolor CT and Targeted Gold Nanoparticles1[END_REF].

Material decomposition is a nonlinear inverse problem [START_REF] Abascal | Nonlinear material decomposition using a regularized iterative scheme based on the Bregman distance[END_REF]. State-of-the art for solving this problem is based on variational approaches [START_REF] Ducros | Regularization of nonlinear decomposition of spectral x-ray projection images[END_REF][START_REF] Abascal | Nonlinear material decomposition using a regularized iterative scheme based on the Bregman distance[END_REF] but the disadvantage of these techniques is that they rely on accurate knowledge of the forward model, including the scanner energy response, which is generally unknown.

Deep learning-based methods are expected to overcome these hurdles. On the one hand, they can learn the inverse problem without a priori assumptions about the scanner energy response. On the other hand, deep convolutional neural networks (CNN) have shown outstanding results in several image processing tasks [START_REF] Lecun | Deep learning[END_REF] and are also promising for learning several inverse problems [START_REF] Hauptmann | Model-Based Learning for Accelerated, Limited-View 3-D Photoacoustic Tomography[END_REF][START_REF] Ravishankar | Image Reconstruction: From Sparsity to Data-adaptive Methods and Machine Learning[END_REF][START_REF] Arridge | Solving inverse problems using data-driven models[END_REF]. Deep This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement N • 701915. It was also performed within the framework of the LabEx PRIMES (ANR-11-LABX-0063) of University de Lyon and under the support of the ANR project SALTO (ANR-17-CE19-0011-01) . The project has also received funding from the European Union's Horizon 2020 research and innovation programme under the grant agreement N 668142. This work was partly funded by France Life Imaging (grant ANR-11-INBS-0006) from the French Investissements d'Avenir. learning has been also proposed for solving the material decomposition in the image domain [START_REF] Touch | A neural network-based method for spectral distortion correction in photon counting x-ray CT[END_REF].

The success of CNNs depends in part on the availability of largescale annotated data sets [START_REF] Krizhevsky | ImageNet classification with deep convolutional neural networks[END_REF]. However, in medical imaging and specifically in emerging imaging modalities such as spectral CT data are very limited. For this, transfer learning appears as relevant technique for medical imaging [START_REF] Shin | Deep Convolutional Neural Networks for Computer-Aided Detection: CNN Architectures, Dataset Characteristics and Transfer Learning[END_REF].

In this work we propose a deep learning approach in the projection domain based on a U-Net architecture [START_REF] Ronneberger | U-Net: Convolutional Networks for Biomedical Image Segmentation[END_REF] and a transfer learning technique that can learn the energy response and compare it to a previously used regularized Gauss-Newton (RGN) method [START_REF] Ducros | Regularization of nonlinear decomposition of spectral x-ray projection images[END_REF] using numerical human thorax phantoms.

METHODS

Forward model

We assume a 2D sensor with P pixels and I energy bins and an object with V voxels made of M materials. We image the object under Θ projections. Let s = (s 1 1,1 , . . . , s θ i,p , . . . , s Θ I,P ) be the measurement vector, where s θ i,p is the data measured in the ith energy bin at the p-th pixel for the θ-th projection, and ρ = (ρ1,1, . . . , ρm,v, . . . , ρM,V ) be the (unknown) mass densities vector, where ρm,v is the mass density for the m-th material at the v-th voxel. We denote

s = G(ρ) (1)
where G represents the forward model that maps ρ onto s. The goal of SCT is to invert [START_REF] Roessl | K-edge imaging in x-ray computed tomography using multi-bin photon counting detectors[END_REF]. The forward model G can be seen as the composition of the linear X-ray transform and a non-linear spectral mixing operator. The X-ray transform X applies to each material independently, i.e.,

am = X(ρm), 1 ≤ m ≤ M (2)
where ρm = (ρm,1, . . . , ρm,v, . . . , ρm,V ) and am = (a 1 m,1 , . . . , a θ m,p , . . . , a Θ m,P ) represent mass density and projected mass density for the m-th material, respectively. Spectral mixing applies to each view independently, i.e.,

s θ = F (a θ ), 1 ≤ θ ≤ Θ (3)
where a θ = (a θ 1,1 , . . . , a θ m,p , . . . , a θ M,P ) and s θ = (s θ 1,1 , . . . , a θ i,p , . . . , a θ I,P ) . In particular, we consider the following non linear mixing [START_REF] Schlomka | Experimental feasibility of multi-energy photon-counting K-edge imaging in pre-clinical computed tomography[END_REF]

s θ i,p = E n0(E)di,p(E) exp - M m=1 a θ m,p τm(E) dE ( 4
)
where E is the energy range of the source, n is the source energy spectrum, di,p is the detector response function at the p-th pixel for the i-th bin, and τm is the m-th material mass attenuation coefficient.

Image reconstruction

We solve the inverse spectral CT image reconstruction problem in two steps: material decomposition and tomographic reconstruction. Material decomposition is the inversion of (3), which is done projection by projection, while tomographic reconstruction is the inversion of (2). In this work, we consider the material decomposition problem with two different approaches, as described below. Assuming parallel geometry, we perform tomographic reconstruction using filtered back-projection algorithm with Ramp filter.

Material decomposition with a variational method

Using a variational framework, material decomposition has been previously formulated as the minimization of the following cost function [START_REF] Ducros | Regularization of nonlinear decomposition of spectral x-ray projection images[END_REF]:

C(a θ ) = 1 2 s θ -F (a θ ) 2 W θ + α m Rm(a θ ), 1 ≤ θ ≤ Θ (5)
where the data fidelity term is a weighted least square functional to approximate a Poisson distribution,

W θ = diag(1/ √ s θ
) is a weighting matrix, α is a regularization parameter, and Rm accounts for specific material regularization. We chose second-order and first-order Tikhonov regularization for soft tissue and bone, respectively, and Huber functional for contrast agent as suggested in [START_REF] Ducros | Regularization of nonlinear decomposition of spectral x-ray projection images[END_REF]. Cost function ( 5) is efficiently minimized by using Gauss Newton algorithm (RGN); MATLAB code for RGN is available from Spectral X-ray image reconstruction SPRAY toolbox [START_REF] Ducros | SPRAY -Spectral X-ray image reconstruction toolbox[END_REF].

Material decomposition by deep learning method

Deep learning has been recently proposed for learning inverse problems [START_REF] K H Jin | Deep Convolutional Neural Network for Inverse Problems in Imaging[END_REF][START_REF] Kang | A deep convolutional neural network using directional wavelets for low-dose X-ray CT reconstruction[END_REF]. In this work, we want to learn

h β : s θ → a θ (6)
where β indicates the parameters of the neural network. We use a U-Net architecture [START_REF] Ronneberger | U-Net: Convolutional Networks for Biomedical Image Segmentation[END_REF] that consists of 14 hidden layers and three levels as shown in figure 1.

Numerical phantom data

Numerical human thorax phantoms are comprised of soft tissue, bone and kidneys marked with gadolinium (Gd) and are built from CT volumes obtained from the kits19 challenge data set [START_REF] Heller | The KiTS19 Challenge Data: 300 Kidney Tumor Cases with Clinical Context, CT Semantic Segmentations, and Surgical Outcomes[END_REF]. The distribution of concentration of gadolinium within the kidneys is generated as homogeneous, where the concentration for each phantom is obtained by sampling from a normal distribution centered in 0.1 g.cm -3 and standard deviation of 0.02 g.cm -3 . Phantoms have 640 × 640 × 100 pixels and after projection each data set consists of 360 projections over a 180 angle span with size 909 × 100 detector pixels. Spectral CT data is simulated using SPRAY toolbox with the modeling parameters (n0(E), τm(E), di(E) (4)) used in [START_REF] Ducros | Regularization of nonlinear decomposition of spectral x-ray projection images[END_REF]. Noisy data is considered by assuming a Poisson distribution for a total number of photons N0 = 10 7 .

To assess the effect of model deviations, we perturb the detector response function considering a pixel-dependent energy shift

di,p(E) = di(E -∆Ei,p), (8)
where ∆Ei,p is the energy shift at pixel p and bin i. This shift is sampled from a normal distribution with a standard deviation of 3 keV. The data simulated with the ideal energy response di(E) (which is constant across detector pixels) are referred to as ideal, while the data simulated with the perturbed energy response di,p(E) are referred to as non-ideal.

For the ideal case, we split phantoms in 40 phantoms for training, 10 for validation and 10 for test, where there are 360 projection images per phantom. We pretain U-Net on ideal data and compare it to RGN. In order to assess the effect of non-ideal data, we assess U-Net on a non-ideal version of the test data. Then, in order to account for deviations from ideality conditions, we applied U-Net with Transfer Learning (U-Net+TL), for which we fine tune the pretrained network (on ideal data) using one phantom for non-ideal data.

RESULTS

Material decomposed projections

Figure 2 shows phantom and decomposed projections. For an ideal detector, RGN leads to excessively smooth images for α = 1.4; smaller values of α led to noisier images. Hence, RGN presents a strong trade-off between noise and regularization. On the contrary, U-Net removes noise while maintaining image quality.

For a non-ideal detector, both methods fail to provide accurate material decomposition (results shown only for U-Net). Assuming perfect data conditions translate into noisy projections. U-Net with transfer learning leads to a significantly reduction in noise.

Tomographic reconstructed images

Figure 3 shows phantom and tomographic reconstructed images. For ideal data, RGN leads to loss of details in soft tissue and bone. U-Net is able to decrease noise while maintaining image quality. For non-ideal data, U-Net presents ring artifacts and large crosstalk between soft tissue and Gd. However, U-Net+TL, which accounts for deviations from the ideal detector assumption, leads to reduction of most artefacts.

Figure 4 shows quantitative results. For ideal data, U-Net leads up to a three fold reduction on MSE on decomposed projections and reconstructions with respect to RGN on the three materials (MSE on the projections not shown). With respect to quantification of Gd, errors were less than 2% with both methods. RGN provides lower mean quantification error but larger variance is present in the histogram of Gd image (the latter not shown).

For non-ideal data, U-Net leads up to 4-fold increase in MSE and 20-fold increase in quantification error. Nevertheless, the proposed U-Net+TL yields 4-fold error reduction.

DISCUSSION

We have proposed a deep learning approach for solving the material decomposition problem that can learn the detector response of a SCT scanner. The approach is based on a U-Net architecture and a transfer learning approach to account for non ideal data.

For ideal data, U-Net provided superior image quality than RGN. This can be explained by the fact that U-Net can learn implicitly the probability distribution of the decomposed materials. In fact, Bayesian deep learning interprets supervised learning as recovering the posterior and it approximates the conditional mean when using the MSE loss [START_REF] Adler | Deep Bayesian Inversion Computational uncertainty quantification for large scale inverse problems[END_REF]. This is particularly relevant for material decomposition in the projection domain, as in this case the choice of prior distributions is less clear. A second disadvantage of variational methods is the high dependence on the regularization parameter.

In terms of quantification of gadolinium, for ideal data, both methods led to less than 2 % quantification error but RGN led to lower error than U-Net. However, inspection of the histogram reveals that RGN presents higher variance. High quantification performance for RGN could be expected as we have considered a low noise scenario and a perfect forward model for RGN. We found that RGN performance decreased for higher noise settings (results not shown).

For non ideal data, methods failed to provide satisfactory material decomposition, with the apparition of large ring artefacts and a 20-fold increase in quantification error. Ring artefacts arise from bias errors in the detector pixels. To mitigate these errors calibration is generally performed in practice. Nevertheless, even after calibration, errors may remain. Here, data discrepancy was modelled as a perturbation of the detector response where the size of the perturbation can be considered as moderate. Lower errors should be encountered in practice after recalibration.

To account for non ideal data, we proposed a transfer learning methodology based on pretraining on numerical phantoms and then fine tuning with one new data set corresponding to non-ideal data. These results show that few data are needed to correct for model deviations. Fine tuning can be understood as recalibrating the model. In practice, a first recalibration could require more data and computational training time yet subsequent recalibrations should need less resources.

There are different approaches to transfer learning. It is common routine during the fine tuning step to freeze some layers by fixing the weights to those values obtained in the pretraining step. However, we preferred to use the pretraining weights as an initial estimate but let them freedom to change [START_REF] Erhan | Why does Unsupervised Pretraining Help Deep Learning?[END_REF]. One reason is that frozen layers are generally the first ones, which capture lower level features, but in our case it could be expected that the first layers are actually those that would overtake the recalibration. In any case, pretraining has been shown to have a regularization (prior) effect, making the solution to lie in a specific region of the solution space, with better generalization properties [START_REF] Erhan | Why does Unsupervised Pretraining Help Deep Learning?[END_REF], so it seems ideal for spectral CT for which data are limited.

This work is subject to few limitations. We have assessed the effect of non-ideal data by perturbing the detector response by a shift in energy. While this model may not be completely exact in practice, it shows the effect of mismodelling the energy response of the scanner and that the proposed approach can learn these type of model deviations. With regard to the fine tuning step, it requires a new 'experimental' data set corresponding to non-ideal data. In practice, this new data set could be an experimental phantom, which properties are known. In this work, we consider a thorax phantom, which is not entirely realistic, but demonstrates the method's capabilities and shows that few data are needed for recalibration. Future work will also consider experimental data.

In conclusion, deep learning has a great potential for SCT as it allows for material decomposition without prior knowledge of the scanner energy response and provides superior image quality than previous variational approaches.

Fig. 1 :

 1 Fig. 1: U-Net architecture. Input is a stack of photon counting projections (I energy bins) and output is a stack of the decomposed material projection images (M materials) for a projection angle.Learning means finding the parameters β that minimizes the following loss function

Fig. 2 :Fig. 3 :

 23 Fig. 2: Left: Phantom and material decomposed projections with RGN and U-Net for ideal data. Right: Phantom and material decomposed projections with U-Net with and without transfer learning for non-ideal data. Projection correspond to an angle view of 0 • .

Fig. 4 :

 4 Fig. 4: MSE of reconstructed images, for soft tissue (S), bone (B) and Gd (G), and Gd quantification error. Left: U-Net (UN) vs RGN (GN). Right: UN vs UN with transfer learning (UNTL).