J P Da Costa

P Weisbecker

B Farbos

J.-M Leyssale

G L Vignoles

C Germain

Investigating carbon materials nanostructure using image orientation statistics

A new characterization method of the lattice fringe images of turbostratic carbons is proposed. This method is based on the computation of their orientation field without explicit detection of fringes. It allows meaningful insights into the material nanostructure and nanotexture at several scales, either qualitatively or quantitatively. The calculation of pairwise spatial statistics of the orientation field at short distance provides measurements of the coherence lengths along any direction, in particular along and orthogonally to the layers. These statistics also allow representing orientation coherence patterns typical of the observed nanostructure. At larger distances, the mean disorientation of the fringes is computed and information about the homogeneity of the sample is obtained. An experimental validation is carried out on various artificial images and an application to the characterization of four bulk turbostratic carbons is provided.

Introduction

Lattice fringe (LF) imaging using High Resolution Transmission Electron Microscopy (HRTEM) has been for more than two decades a choice technique for investigating and quantifying the nanostructure (i.e. the spatial arrangements of graphene layers) of carbonaceous materials such as soot [START_REF] Fernandez-Alós | Soot and char molecular representations generated directly from HRTEM lattice fringe images using Fringe3D[END_REF], activated carbons [START_REF] Pré | A new approach to characterize the nanostructure of activated carbons from mathematical morphology applied to high resolution transmission electron microscopy images[END_REF], chars [START_REF] Shim | A methodology for analysis of 002 lattice fringe images and its application to combustion-derived carbons[END_REF][START_REF] Mathews | The utility of coal molecular models[END_REF][START_REF] Vander Wal | Analysis of HRTEM images for carbon nanostructure quantification[END_REF], pyrocarbons [START_REF] Reznik | On the terminology for pyrolytic carbon[END_REF][START_REF] Bourrat | Regenerative laminar pyrocarbon[END_REF][START_REF] Leyssale | Structural features of pyrocarbon atomistic models constructed from transmission electron microscopy images[END_REF][START_REF] Farbos | Nanoscale structure and texture of highly anisotropic pyrocarbons revisited with transmission electron microscopy, image processing, neutron diffraction and atomistic modelling[END_REF], or nuclear graphites [START_REF] Ammar | Characterization of graphite implanted with chlorine ions using combined Raman microspectrometry and transmission electron microscopy on thin sections prepared by focused ion beam[END_REF][START_REF] Gotoh | High resolution electron microscopy of graphite defect structures after keV hydrogen ion bombardment[END_REF].

Though meaningful and informative, the qualitative interpretation of LF images is subjective and time consuming. The need for quantitative analysis tools has rapidly lead to develop automated image analysis tools, as reported e.g. in refs. [START_REF] Shim | A methodology for analysis of 002 lattice fringe images and its application to combustion-derived carbons[END_REF][START_REF] Sharma | A new quantitative approach for microstructural analysis of coal char using HRTEM images[END_REF][START_REF] Rouzaud | Quantitative high-resolution transmission electron microscopy: a promising tool for carbon materials characterization[END_REF][START_REF] Yehliu | Development of an HRTEM image analysis method to quantify carbon nanostructure[END_REF].

Whatever the type of carbon materials under study, most of published LF image analysis procedures are based on four main algorithmic steps, namely (i) filtering of the digitized image, (ii) fringe extraction, (iii) morphological analysis of individual fringes, (iv) characterization of spatial arrangements of fringes. Essentially, image filtering (i) aims at removing high frequency artifacts (mainly electronic noise) and low frequency variations (Bragg artifacts and sample thickness variations) that are not meaningful to describe material nanostructure and may hinder image analysis. Filtering is usually conducted in the Fourier domain by applying either ideal square filters [START_REF] Pré | A new approach to characterize the nanostructure of activated carbons from mathematical morphology applied to high resolution transmission electron microscopy images[END_REF][START_REF] Shim | A methodology for analysis of 002 lattice fringe images and its application to combustion-derived carbons[END_REF][START_REF] Rouzaud | Quantitative high-resolution transmission electron microscopy: a promising tool for carbon materials characterization[END_REF] or smooth Gaussian filters [START_REF] Tóth | Quantitative differentiation of poorly ordered soot nanostructures: a semi-empirical approach[END_REF]. As mentioned by Toth and coworkers, the latter prevents harmonic ringing caused by the sharp changes in the transfer function and subsequent false fringe detections. In the case of more anisotropic materials such as rough or regenerated laminar pyrocarbons for instance [START_REF] Leyssale | Structural features of pyrocarbon atomistic models constructed from transmission electron microscopy images[END_REF], a combination of radial and angular filtering has also appeared to be a relevant technique to remove irrelevant details in LF images.

Typically, fringe detection (ii) starts with image binarization in order to discriminate fringe objects from background. Binarization is obtained by thresholding using either a fixed [START_REF] Shim | A methodology for analysis of 002 lattice fringe images and its application to combustion-derived carbons[END_REF][START_REF] Rouzaud | Quantitative high-resolution transmission electron microscopy: a promising tool for carbon materials characterization[END_REF] or an adaptive [START_REF] Tóth | Quantitative differentiation of poorly ordered soot nanostructures: a semi-empirical approach[END_REF] threshold. Rouzaud and Clinard [START_REF] Rouzaud | Quantitative high-resolution transmission electron microscopy: a promising tool for carbon materials characterization[END_REF] and later Pré et al. [START_REF] Pré | A new approach to characterize the nanostructure of activated carbons from mathematical morphology applied to high resolution transmission electron microscopy images[END_REF] suggested the prior use of a top-hat transform to ease the thresholding process. Binarization is usually followed by fringe skeletonization and various post-processing operations in order to separate, reconnect or remove fringe chunks [START_REF] Pré | A new approach to characterize the nanostructure of activated carbons from mathematical morphology applied to high resolution transmission electron microscopy images[END_REF][START_REF] Shim | A methodology for analysis of 002 lattice fringe images and its application to combustion-derived carbons[END_REF][START_REF] Sharma | A new quantitative approach for microstructural analysis of coal char using HRTEM images[END_REF][START_REF] Yehliu | Development of an HRTEM image analysis method to quantify carbon nanostructure[END_REF]. Another possible approach for fringe detection relies in the use of a sub-pixel level curve tracking algorithm initially proposed in [START_REF] Costa | Level Curve Tracking Algorithm for Textural Feature Extraction[END_REF] and later used in [START_REF] Costa | A Curvilinear Approach for Textural Feature Extraction: Application to the Characterization of Composite Material Images[END_REF][START_REF] Raynal | Multi-scale quantitative analysis of carbon structure and texture: III. Lattice fringe imaging analysis[END_REF]. After image background variations removal, level curves are tracked on either side of the fringes around the mean intensity level until the fringe ends (zones with poor orientation coherence) are reached.

Morphological fringe analysis (iii) usually comprises statistics of individual fringe lengths and tortuosities [START_REF] Pré | A new approach to characterize the nanostructure of activated carbons from mathematical morphology applied to high resolution transmission electron microscopy images[END_REF][START_REF] Shim | A methodology for analysis of 002 lattice fringe images and its application to combustion-derived carbons[END_REF][START_REF] Yehliu | Development of an HRTEM image analysis method to quantify carbon nanostructure[END_REF][START_REF] Raynal | Multi-scale quantitative analysis of carbon structure and texture: III. Lattice fringe imaging analysis[END_REF], curvatures and orientations [START_REF] Pré | A new approach to characterize the nanostructure of activated carbons from mathematical morphology applied to high resolution transmission electron microscopy images[END_REF][START_REF] Tóth | Quantitative differentiation of poorly ordered soot nanostructures: a semi-empirical approach[END_REF]. Shim et al. also proposed a measure of the apparent crystalline fraction [START_REF] Shim | A methodology for analysis of 002 lattice fringe images and its application to combustion-derived carbons[END_REF].

Besides individual fringe features statistics, various authors addressed the description of mutual arrangements of fringes (iv) by measuring interlayer spacing histograms [START_REF] Pré | A new approach to characterize the nanostructure of activated carbons from mathematical morphology applied to high resolution transmission electron microscopy images[END_REF][START_REF] Rouzaud | Quantitative high-resolution transmission electron microscopy: a promising tool for carbon materials characterization[END_REF][START_REF] Yehliu | Development of an HRTEM image analysis method to quantify carbon nanostructure[END_REF].

Toth et al. [START_REF] Tóth | Quantitative differentiation of poorly ordered soot nanostructures: a semi-empirical approach[END_REF] also proposed to measure the interlayer spacing deviation parameter or the junction parameter (number of branch points per unit length). The identification of stacks of parallel fringesassimilated to Basic Structural Units [START_REF] Oberlin | High-resolution TEM studies of carbonization and graphitization[END_REF][START_REF] Rouzaud | Structure, microtexture, and optical properties of anthracene and saccharose-based carbons[END_REF] has also been addressed by several authors [START_REF] Pré | A new approach to characterize the nanostructure of activated carbons from mathematical morphology applied to high resolution transmission electron microscopy images[END_REF][START_REF] Rouzaud | Quantitative high-resolution transmission electron microscopy: a promising tool for carbon materials characterization[END_REF][START_REF] Yehliu | Development of an HRTEM image analysis method to quantify carbon nanostructure[END_REF][START_REF] Tóth | Quantitative differentiation of poorly ordered soot nanostructures: a semi-empirical approach[END_REF]. They established statistics for the stacking number N and the coherence length L a . Other valuable features that can be extracted from LF images are the 2D nematic and polar order parameters [START_REF] Shim | A methodology for analysis of 002 lattice fringe images and its application to combustion-derived carbons[END_REF]. Computed from the orientations of the fringes within a given spatial window around any given location of the image, these order parameters quantify the organization of fringe orientations inside the window. In brief, the nematic order parameter gets high values where fringes are all nearly parallel while the polar parameter is high in the case of concentric fringes around the window center. Such parameters have later been used as saliency features to detect soot spherules in large images and drive an automated image analysis process [START_REF] Tóth | Automated analysis of heterogeneous carbon nanostructures by high-resolution electron microscopy and online image processing[END_REF].

Though yielding lots of relevant information about many material nanostructures, the bottleneck of this kind of approaches is the need for a fringe extraction algorithm that usually consists of numerous complex steps with various parameters, sometimes subjective and most of the time application specific. In recent works, some authors argue that valuable information can also be obtained without explicit fringe extraction. For instance, Toth et al. [START_REF] Tóth | A novel framework for the quantitative analysis of high resolution transmission electron micrographs of soot I. Improved measurement of interlayer spacing[END_REF] proposed an original framework to compute statistics of interlayer spacing and fringe orientation based on the direct application of a Gabor filter bank on an LF image. In [START_REF] Tóth | A novel framework for the quantitative analysis of high resolution transmission electron micrographs of soot II. Robust multiscale nanostructure quantification[END_REF], the same authors presented an alternative approach for the computation of multiscale nematic and polar order parameters. In this approach, the orientation statistics were computed from the image orientation field estimated using steerable filters [START_REF] Freeman | The design and use of steerable filters[END_REF] without explicit identification of fringes.

In former works, Germain et al. [START_REF] Germain | Multiscale estimation of vector field anisotropy -application to texture characterization[END_REF] also proposed a general framework for the multiscale evaluation of texture anisotropy using feature I so which was used to discriminate between various raw and heat treated pyrocarbons.

In this paper, we address the characterization of spatial arrangements of fringes in turbostratic carbons through LF images and their orientation fields. The latter are computed by standard computationally efficient estimation techniques without explicit detection of fringes. We show that the visualization of the local orientation map and of the orientation coherences allows a more convenient representation of the nanotexture and of the defects. In addition, spatial statistics of the orientation field are used to gather information on the texture at several scales. While previous approaches based on the multiscale nematic order [START_REF] Shim | A methodology for analysis of 002 lattice fringe images and its application to combustion-derived carbons[END_REF][START_REF] Tóth | A novel framework for the quantitative analysis of high resolution transmission electron micrographs of soot II. Robust multiscale nanostructure quantification[END_REF] or on multiscale anisotropy [START_REF] Germain | Multiscale estimation of vector field anisotropy -application to texture characterization[END_REF] provide global estimations of the degree of anisotropy at one or several scales, the specificity of our work is to consider pairwise statistics, i.e. statistics of pairs of pixels. Using ad-hoc graphical representations, our approach allows monitoring how fringe orientation coherence vanishes between two sites according to distance and orientation, thus providing access to the grain shapes. The relevance of our approach is first assessed on artificial images showing grain tessellations before studying real material images. In particular, for highly textured materials as turbostratic carbons, we obtain coherence lengths at short distance, parallel and perpendicular to the layers, which are related to other coherence measures provided by X-ray diffraction (XRD). At larger distances, the mean angular deviation of the fringes is computed and information about the homogeneity of the sample is obtained.

The paper is organized as follows. The second section provides details about material samples, characterization techniques and image processing algorithms including orientation retrieval and spatial statistic computation. The third section deals with experimental validation on various artificial images and shows an application to the characterization of four types of bulk carbons. Finally, we will draw some conclusions and mention a few prospects for future work.

Materials and methods

Material samples

Four materials were chosen to evaluate this method: three pyrocarbons (PyCs)namely the Rough Laminar PyC (RL), the Regenerated Laminar PyC (ReL) and the Smooth Laminar PyC (SL)and a PAN-based carbon fiber (C Fiber ). All materials are made of graphene layers stacked together with a turbostratic disorder between the layers. Differences between the materials lay in their texture (i.e. the preferential orientation of the layers), the curvature of their layers, and the shape and size distribution of their BSUs (Basic Structural Units).

The RL PyC [26] and the ReL PyC [START_REF] Bourrat | Regenerative laminar pyrocarbon[END_REF] are both high textured PyCs [START_REF] Reznik | On the terminology for pyrolytic carbon[END_REF], the latter containing a higher amount of structural defects. The smooth laminar PyC (SL) is medium textured and more heterogeneous as it is made of a mixture of lamellae and pores of various sizes [START_REF] Goma | Characterization of low temperature pyrocarbons obtained by densification of porous substrates[END_REF].

C Fiber exhibits a distinct texture, isotropic when observed along the fiber axis and rather homogeneous with strongly folded layers, corresponding to the axial orientational order [START_REF] Shim | A methodology for analysis of 002 lattice fringe images and its application to combustion-derived carbons[END_REF].

ReL, RL and SL PyCs were prepared by isothermal isobaric Chemical Vapor Infiltration (CVI) from respectively pure propane (ReL) and a methane/propane mixture (RL) precursors.

The RL PyC was obtained by deposition on the surface of a fused quartz tube filled at both ends by a carbon preform and placed in the CVI reactor [START_REF] Weisbecker | Microstructure of pyrocarbons from pair distribution function analysis using neutron diffraction[END_REF] while the SL PyC was obtained in the same conditions, but with the fused quartz tube left open. For the ReL PyC, a silica woven fabric was infiltrated using pure propane with the following processing parameters: total pressure p = 5 kPa, temperature T = 1050°C, residence time t r = 3 seconds and 10 hours infiltration duration [START_REF] Weisbecker | A Pair Distribution Function approach to pyrocarbon structure using neutron diffraction[END_REF]. The carbon fiber is a commercial ex-PAN carbon fiber [START_REF] Guigon | Microtexture and Structure of Some High-Modulus, PAN-Based Carbon Fibres[END_REF].

Material characterization 2.2.1 HRTEM imaging

Before being studied by TEM, bulk samples were embedded in epoxy glue, mechanically polished down to 100 m using a 15 m diamond plate and then ion thinned in order to achieve electron transparency using a JEOL Ion Slicer (EM-09100IS) [START_REF] Weisbecker | Thin-film preparation of C/C composites and CMC using the broad argon ion beam method[END_REF]. The initial incidence angle was set to 0.5° with an energy of 5.5 keV for the Ar + ion beam; then, the angle was increased to 2.5° and the energy progressively reduced to 4 keV and eventually 2 keV in order to reduce the irradiation damage to the carbon structure.

High-resolution transmission electron microscopy (HRTEM) images were acquired on a Philips CM30ST TEM (300 kV -LaB 6 ). All images were obtained near Scherzer defocus (-57nm) with an objective aperture of 4 nm -1 allowing to select only the 002 lattice fringes.

Near Scherzer defocus the contrast transfer function (CTF) of the transmission electron microscope (TEM) is nearly flat for an extended region of spatial frequencies (roughly from 2 to 4.2 nm -1 ), including the frequency domain corresponding to the 002 lattice fringes (2.9 nm -1 ). In these conditions the atoms appear in black in the LF image which can be seen, as a first approximation, as an inverted image of the projected potential (projection along the TEM optical axis through the thin foil thickness). It is well known that graphitic structures are prone to amorphization under a 300kV electron beam [START_REF] Tanabe | Radiation damage of graphite -Degradation of material parameters and defect structures[END_REF], hence the acquisition time was reduced to 10 seconds for each zone in order to reduce the irradiation dose as much as possible.

X-ray diffraction

RL and ReL PyCs were analyzed by powder diffraction using a Bruker D8 Advance diffractometer (K  Cu) in order to determine the coherence lengths along the c-axis L c and in the plane L a . Powders were filled into a zero background sample holder and the data were acquired on a 2 range of 10-74°. The whole patterns were fitted using the CarbonXS software [START_REF] Shi | Structure-refinement program for disordered carbons[END_REF] which accounts for the strong overlap between the 10 band and the 004 peak.

The SRM1976 (corundum) calibration standard was used to check the alignment of the diffractometer and to determine the instrumental resolution which is smaller than 0.07 °2 while carbon peaks are larger than 1.6 °2 . Thus the instrumental resolution was not taken into account as it is negligible comparing to the pyrocarbons peaks widths.

Image processing

LF images belong to the category of what is called textured images. Texture is a keystone of many image analysis tasks including classification, segmentation or characterization. In order to extract texture descriptors, statistical approaches based on second and higher order statistics have been widely reported in literature. For instance, methods based on the autocorrelation and on co-occurrence matrices [START_REF] Haralick | Statistical and structural approaches to textures[END_REF] or on gray level difference histograms [START_REF] Unser | Sum and difference histograms for texture classification[END_REF][START_REF] Chetverikov | Texture analysis using feature-based interaction maps[END_REF] are simple and efficient to extract texture descriptors. These methods rely only on the distribution of pixel intensity values. Several authors [START_REF] Davis | Texture analysis using generalized co-occurrence matrices[END_REF][START_REF] Kovalev | Multidimensional co-occurence matrices for object recognition and matching[END_REF][START_REF] Kovalev | Using orientation tokens for object recognition[END_REF] pointed out that methods based on gray level statistics are not suited to the description of structured textures. Indeed, textures arising from an arrangement of local patterns cannot be described only in terms of gray level variations.

Descriptors related to the spatial arrangement of texture patterns are needed. Pattern orientations have been recognized by several authors as important features for machine vision systems [START_REF] Germain | Multiscale estimation of vector field anisotropy -application to texture characterization[END_REF][START_REF] Chetverikov | Texture analysis using feature-based interaction maps[END_REF][START_REF] Davis | Texture analysis using generalized co-occurrence matrices[END_REF][START_REF] Kovalev | Multidimensional co-occurence matrices for object recognition and matching[END_REF][START_REF] Kovalev | Using orientation tokens for object recognition[END_REF][START_REF] Kass | Analysing oriented patterns[END_REF][START_REF] De Luis-Garcia | Texture and color segmentation based on the combined use of the structure tensor and the image components[END_REF]. Although the extraction and the regularization of pixel orientations have been widely addressed, e.g. in [START_REF] Kass | Analysing oriented patterns[END_REF][START_REF] Bigün | Multidimensional orientation estimation with applications to texture analysis and optical flow[END_REF][START_REF] Pouliquen | A new adaptive framework for unbiased orientation estimation[END_REF][START_REF] Michelet | Estimating local multiple orientations[END_REF], very few studies have focused on the use of spatial statistics of the orientation field. In their pioneering works, Davis et al. [START_REF] Davis | Texture analysis using generalized co-occurrence matrices[END_REF] presented a generalization of co-occurrence matrices accounting for pattern orientations. They defined the generalized co-occurrence matrix as a multidimensional array where each element is a count of the number of occurrences of a given spatial configuration of pixels. This configuration includes the relative position of the pixels and any kind of features characterizing the pixels.

The generalized co-occurrences were applied to the magnitudes and orientations of the gradient vectors, computed on object boundaries. Later, Kovalev et al. [START_REF] Kovalev | Multidimensional co-occurence matrices for object recognition and matching[END_REF][START_REF] Kovalev | Using orientation tokens for object recognition[END_REF] provided implementations of multidimensional co-occurrence matrices. They proposed either to combine gray level and orientation statistics [START_REF] Kovalev | Multidimensional co-occurence matrices for object recognition and matching[END_REF] or to base the co-occurrences on the orientation field only [START_REF] Kovalev | Using orientation tokens for object recognition[END_REF]. Their approach, applied to a small set of spatial relations, proved to be efficient for the discrimination of specific textures with local anisotropies.

In the following we propose an original framework for the second order statistical analysis of the orientation field. Based on a rotationally invariant implementation of orientation spatial differences, it specifically addresses the case of anisotropic textures such as those appearing in HRTEM LF images.

Image filtering

HRTEM image filtering is carried out according to the radial-angular filter used, though not detailed, in [START_REF] Leyssale | Structural features of pyrocarbon atomistic models constructed from transmission electron microscopy images[END_REF][START_REF] Farbos | Nanoscale structure and texture of highly anisotropic pyrocarbons revisited with transmission electron microscopy, image processing, neutron diffraction and atomistic modelling[END_REF][START_REF] Leyssale | An image guided atomistic reconstruction of pyrolytic carbons[END_REF]. The filter performs in the Fourier domain and is defined by its band-pass transfer function:

H(𝜌, 𝜃) = H radial (𝜌)H angular (𝜃) with H radial (𝜌) = 1 𝜎 𝜌 √2𝜋 𝑒 - (𝜌-𝜌 0 ) 2 2𝜎 𝜌 2 and H angular (𝜃) = 1 𝜎 𝜃 √2𝜋 𝑒 - (𝜃-𝜃 0 ) 2 2𝜎 𝜃 2
(𝜌, 𝜃) are the polar frequency coordinates in the Fourier domain. 𝜌 0 = 1 𝑑 002 ≈ 2.9 nm -1 corresponds with the interlayer spacing. 𝜎 𝜌 = 𝛼𝜌 0 is defined from 𝜌 0 with 0 < 𝛼 < 1. H radial performs a Gaussian filtering around 𝜌 = 𝜌 0 , reducing high and low frequency artifacts while preserving 𝑑 002 fringes whatever their orientation. The second term H angular aims at preserving patterns which are oriented around a reference orientation 𝜃 0 that fits the preferential orientation in the image. 𝜃 0 is found by looking for the largest magnitude peak in the Fourier transform modulus. Filtering is done using a Gaussian function centered on 𝜃 0 with a standard deviation of 𝜎 𝜃 . The typical values 𝛼 = 0.5 and 𝜎 𝜃 = 50° have proved to be a good compromise between artifact removal and structure preservation, (see e.g. [START_REF] Leyssale | Structural features of pyrocarbon atomistic models constructed from transmission electron microscopy images[END_REF][START_REF] Farbos | Nanoscale structure and texture of highly anisotropic pyrocarbons revisited with transmission electron microscopy, image processing, neutron diffraction and atomistic modelling[END_REF]), for anisotropic PyCs such as the RL and ReL PyCs. For poorly anisotropic materials such as the SL PyC and C Fiber , angular filtering is deactivated in order to prevent the apparition of unwanted artificial structures in the filtered image. Naturally, though guided by experience, the determination of 𝛼 and 𝜎 𝜃 remains arbitrary and is subject to the appreciation of the experimenter. Nevertheless, when kept within acceptable ranges leading to eye-friendly images, the choice of values is not critical provided that the values are kept constant when comparing various material images together. Indeed, variations in the filter strength or in the filter nature do not hinder the calculation of orientations statistics contrary to fringe detection algorithms that directly depend on image grey level variations. A comparison of various filtering approaches is provided in the supplementary material, showing their effect on the structural characteristics extracted from the orientation field. It is shown that in the absence of filteringor equivalently in the presence of noisethe orientation maps can still be estimated and their statistical analysis still provides information about orientation spatial coherence.

However, the presence of high frequency noise may lead to slightly underestimate the material anisotropy, which can be easily avoided by filtering.

Orientations: definition, estimation and representation

The second step of the LF image processing process is the estimation of the orientation field.

Various definitions of orientation are available in image processing literature. For instance, the linearly symmetric model [START_REF] Bigün | Multidimensional orientation estimation with applications to texture analysis and optical flow[END_REF], the i1D case [START_REF] Krieger | Nonlinear image operators for the evaluation of local intrinsic dimensionality[END_REF] or the mono-directional model [START_REF] Michelet | Estimating local multiple orientations[END_REF] provide definitions we can refer to. They all consider that a unique orientation is available at any location, i.e. it is the orientation of the underlying texture pattern.

In the case of locally unique orientations, a standard approach to extract local orientations in an image is the local structure tensor (LST) [START_REF] Kass | Analysing oriented patterns[END_REF], based itself on the image intensity gradient.

Let I be the image magnitude and Ω  ℕ² the image pixel grid. The image intensity gradient at location 𝐮 ∈ 𝛺 is defined by:

𝛁𝐼(𝐮) = [ 𝐼 𝑥 (𝐮) 𝐼 𝑦 (𝐮) ] = [ 𝜕𝐼(𝐮) 𝜕𝑥 𝜕𝐼(𝐮) 𝜕𝑦 ]
where 𝐼 𝑥 and 𝐼 𝑦 are the partial derivatives of the image. Gradient estimation is usually done using convolution masks in the spatial domain such as local derivative operators (e.g. [START_REF] Pouliquen | A new adaptive framework for unbiased orientation estimation[END_REF]) or scalable gradient filters as Deriche's filter or Gaussian derivatives.

The LST at location 𝐮 is the local covariance of the gradient field:

𝐓(𝐮) = [ T xx (𝐮) T xy (𝐮) T xy (𝐮) T yy (𝐮) ] = 〈𝛁(𝐮)𝛁 † (𝐮)〉 𝑊 𝐮 = [ 〈𝐼 𝑥 ²(𝐮)〉 〈𝐼 𝑥 (𝐮)𝐼 𝑦 (𝐮)〉 〈𝐼 𝑥 (𝐮)𝐼 𝑦 (𝐮)〉 〈𝐼 𝑦 ²(𝐮)〉 ]
where † denotes matrix transpose and The size of the averaging window 𝑊 𝐮 for the LST computation must be tuned according to the expected scale of the orientation estimation. In other words it must fit the size of the patterns of interest. In particular, in the case of lattice fringe images, the size of the averaging window must be tuned to the width of a fringe, the d 002 distance. Regarding the gradient estimation prior to the LST computation, any operator based on convolution filters may be chosen provided that their spatial extent is small compared to the size of the LST averaging window. For instance, 3 × 3 convolutive filters such as Sobel's, Prewitt's or GOP3 [START_REF] Pouliquen | A new adaptive framework for unbiased orientation estimation[END_REF] are eligible.

Many other orientation estimation methods exist including those based on the use of filter banks such as the IRON operators [START_REF] Michelet | Estimating local multiple orientations[END_REF] or the steerable filters [START_REF] Freeman | The design and use of steerable filters[END_REF] used in [START_REF] Tóth | A novel framework for the quantitative analysis of high resolution transmission electron micrographs of soot II. Robust multiscale nanostructure quantification[END_REF] for anisotropy analysis on soot images. The main advantage of filter bank approaches is their ability to detect and discriminate between several orientations simultaneously present on a site [START_REF] Michelet | Estimating local multiple orientations[END_REF] as it may occur in LF images of soot, coal chars or activated carbons images for instance. In most parts of turbostratic carbon images, fringes can be considered as elongated patterns with a unique orientation and without pattern superimpositions. In this case, the use of the structure tensor, much more efficient computationally speaking and far simpler to implement, is more relevant.

Let us note that pattern orientations can be handled as axial data, i.e. 𝜋-periodic circular data [START_REF] Mardia | Directional Statistics[END_REF]. Indeed two opposite angles, 𝜃 0 and 𝜃 1 = 𝜃 0 + 𝜋, relate to the same pattern orientation.

Consequently, orientations 𝜃(𝐮) retrieved from the LST can be restricted to an interval of width 𝜋 such as [0, 𝜋[.

The couple composed of 𝜃(𝐮) and 𝜂(𝐮)(also noted 𝜃 𝐮 and 𝜂 𝐮 ), can be represented as a vector field 𝐕 𝐮 = (𝜂 𝐮 cos 𝜃 𝐮 , 𝜂 𝐮 sin 𝜃 𝐮 ) † from which it is possible to extract spatial pairwise statistics as we will see below.

Spatial statistics of orientations

An intuitive way to measure the spatial dependence of orientations is to compute spatial orientation differences. In order to deal with the axial nature of orientation data, some of us [START_REF] Costa | Orientation difference statistics for texture description[END_REF] proposed to measure the difference between two orientations (𝜃 1 , 𝜃 This feature allows measuring orientation dependencies when moving in a specific direction regarding the local structure orientation as shown in figure 1. for instance the feature-based interaction maps which are intensity coded representations of a spatial feature as a function of the spacing vector [START_REF] Chetverikov | Texture anisotropy, symmetry, regularity: recovering structure and orientation from interaction maps[END_REF].

Regarding computational aspects, the approach proceeds in three major steps; (i) image filtering, (ii) calculation of the orientation and coherence maps and (iii) calculation of the pairwise statistics. For an image of size 𝑁 = 𝑁 𝑥 × 𝑁 𝑦 , step (i) is performed very efficiently with a computational complexity of 𝑂(𝑁log𝑁) in the Fourier domain. The complexity of step (ii) is 𝑂(𝑁 × 𝑊 𝑥 × 𝑊 𝑦 ) , where 𝑊 𝑥 and 𝑊 𝑦 are the dimensions of the tensor averaging window. It is linear in 𝑁, allowing very low computational needs. Steps (i) and (ii) can thus be performed in a few seconds in any personal computer. In comparison, step (iii) is more time consuming as it involve a computational complexity of 𝑂(𝑁 × 𝑁  × 𝑁 𝜌 ) where 𝑁  and 𝑁 𝜌 are the number of angles and distances to be investigated in the pairwise statistics.

Typically, for an image of size 76nm × 76nm, acquired at a resolution 𝑟 = 0.037 𝑛𝑚/𝑝𝑖𝑥𝑒𝑙, the computation of 𝑅𝐼𝑀𝑂𝐷(𝜌, ) for 𝑁  = 2 angles ( = 0° and 90°) and 𝑁 𝜌 = 67 distances (between 0 and 5 𝑛𝑚, every 0.074 𝑛𝑚) takes around 10 minutes on a 2.5 GHz processor, without parallelization.

Results and discussion

Validation of the method on artificial images

To illustrate and evaluate which kind of quantitative information can be provided by orientation statistics, artificial lattice fringe images were generated. The latter are Voronoi tessellations filled with 002 fringes spaced by 0.34 nm. A tessellation is obtained by first generating a given number of seeds randomly in the image and, second, by assigning each pixel to the closest seed, in the sense of the Euclidean distance. The shapes of the Voronoi domains are not explicitly defined but are the result of the assignment process. Frontiers between neighboring domains are linear segments. Note that the shapes of the domains can be easily modified by weighting horizontal and vertical coordinates differently in the distance definition. For instance, if the horizontal coordinate differences are weighted twice as much as the vertical ones, the tessellation process results in rather vertical domains with a shape anisotropy factor of two (here a Y/X growth ratio of two).

Domain orientations are realizations either of (i) a wrapped normal random variable [START_REF] Mardia | Directional Statistics[END_REF] with zero mean (horizontal fringes) and with a given standard deviation, or Plots of feature 𝑅𝐼𝑀𝑂𝐷(𝜌, ) for both images, given at fig. 3, show the orientation difference as a function of distance for  = 0° (along the fringes) and  = 90° (orthogonally to the fringes). Both diagrams have a similar shape: a quick loss of orientation at short distance and a plateau reached at roughly 2.5 nm. However the maximum value of orientation difference is around 12° in the textured sample and 45° in the isotropic sample. The loss of orientation at short distance is due to the interaction between neighboring grains and the plateau corresponds to the mean disorientation between grains. This plateau is reached for distances larger than the largest domains, thus giving domain size information. In order to obtain coherence lengths that can be accurately determined we have selected two parameters corresponding to the distance at which 90% of the value of the plateau is reached for the diagrams at 0° and 90°. These parameters 𝐿 𝑎𝑀𝑂𝐷 and 𝐿 𝑐𝑀𝑂𝐷 are good estimates of the domain sizes as it is shown in table 1. It is also shown that the match between the domain size and 𝐿 𝑐𝑀𝑂𝐷 and 𝐿 𝑎𝑀𝑂𝐷 is good for domain sizes ranging from 0.5 to 3 nm. The value of the plateau, 𝛽 𝑀𝑂𝐷 is a quantification of the mean twist disorientation and provides information similar to the orientation angle (OA°) obtained by electron diffraction but at a lower scale. From table 1 it can be see that 𝛽 𝑀𝑂𝐷 increases when the domain size decreases (standard deviation of orientation =9° for all images); this is due to the higher amount of grain boundaries whose orientation is usually more difficult to determine. The ratio 𝐿 𝑐𝑀𝑂𝐷 /𝐿 𝑎𝑀𝑂𝐷 appears to be slightly larger than unity, suggesting a better orientation coherence orthogonally to the fringes. It can be attributed to the phase shift of the fringes of two neighboring grains having similar orientation; this leads to sharp transitions at the grain boundaries in the fringe direction while transitions are smoother in the direction perpendicular to the fringes. 𝑛 is the number of domains in the image and 𝐷 is the mean domain size.

Experimental images usually exhibit domains with an anisotropic shape. In order to study the efficiency of our method with such a microstructure, we have generated artificial images with columnar grains and with fringes preferentially oriented perpendicular to the columns. In order to demonstrate the ability of the method to gather information on spatial coherence of orientation at both short and long distances we have generated an image corresponding to a loose stacking of concentric circles (concentric cylinders or spheres seen in cross section, that can be considered as rough models of ex-PAN carbon fiber or soot structures), the external circles being 6 nm in diameter. This image as well as the corresponding local orientation map and the mean orientation difference (𝑅𝐼𝑀𝑂𝐷(𝜌, )) are shown in figure 5. The 𝑅𝐼𝑀𝑂𝐷(𝜌, )

plot shows two main features. First at short distances a quick loss of orientation is observed for  = 0° while the slope is smaller for  = 90°, which is due to the perfect stacking of the layer in the radial direction while along the circumference the loss of orientation is due to the curvature of the layers, a maximum disorientation is obtained at 2.6 nm which corresponds to the mean radius of curvature. Second the 𝑅𝐼𝑀𝑂𝐷(𝜌, ) plots shows an oscillating behavior at larger distances, this is due to the stacking of the circles, the first minimum is reached at approximately 6 nm which is the mean distance between two circles.

Application to lattice fringe images

In this section, we illustrate how the spatial statistics of orientations can be meaningfully exploited for the characterization of material images. The materials we are interested in are four dense carbons observed by high resolution transmission electron microscopy (HRTEM).

Four images are represented in figure 6. Resolution is 0.037 nm/pixel. Image size is 2048×2048 pixels i.e. about 76×76 nm 2 (for clarity reasons, only zooms of 38 nm and 9 nm are shown). Depending on the material, the spatial structure of the fringes is more or less regular. For instance, the regenerative laminar pyrocarbon (ReL) and the rough laminar pyrocarbon (RL) are very anisotropic both at low and high scales. Beside their overall anisotropy, a close visual analysis reveals a local organization in stacks of near-parallel fringes. These stacks appear to be larger in RL than in ReL. On the contrary, ex-PAN carbon fiber C Fiber and Smooth Laminar PyC (SL) are anisotropic at a very low scale and show much more disorder at higher scales. C Fiber fringes are organized in compact wrapped patterns while orientations of fringes in SL images show non stationarities at high scale as indicated in the orientation fields of figure 7. Orientation confidences are shown figure 8 for the four materials, these maps enlighten areas of defects (low confidence) when fringes are too disrupted to allow a good determination of the orientation: dislocations, fringe overlapping, pores, and amorphous areas. The two high textured pyC images exhibit few defects which mainly consist of line defects in the images, few nanometers long, oriented mostly perpendicular to the fringes and located at domains boundaries. The amount of these defective areas is clearly larger for ReL than for RL. SL orientation coherence maps show a high amount of defective areas which are partly due to some fringe overlapping. Indeed, the materials being more isotropic, fringes with different orientations can be found through the thickness of the TEM sample. Note however that fringe overlap has only little effect on the RIMOD statistic as these areas exhibit low confidence values (dark areas in the orientation maps, figure 8) which are used to weight the data.

Besides, numerous pores are present inside the sample due to the high curvature of the fringes. These pores are also present in C Fiber material [START_REF] Guigon | Microtexture and Structure of Some High-Modulus, PAN-Based Carbon Fibres[END_REF] as well as some fringe overlapping leading to a relatively high amount of defects in the orientation confidence maps. These pores and the disordered areas appear as areas of weak contrast in the HRTEM images. Orientation fields were obtained on the four images after frequential filtering. Bandpass radial filtering (𝜌 0 = 2.9 nm -1 , 𝛼 = 1) was used for the SL PyC and the CFiber; both radial and angular (𝛼 = 0.5, 𝜎 𝜃 = 50°) filtering for the two high textured PyCs (RL and ReL) in order to reduce various artifacts such as Bragg fringes or Moirés. It must be pointed out that the image filtering is not mandatory to carry out the analysis and a comparative study of various filters is available as supplementary material. It shows that filtering for high textured materials improves the accuracy of the determination of the orientation at the domain boundaries.

When looking at the four orientation fields, it is rather easy to understand that the medium and high scale layout of such locally parallel textures is much easier to analyze by studying spatial orientation statistics. In particular, the use of the rotation invariant feature 𝑅𝐼𝑀𝑂𝐷(𝜌, ) is required here not only because sample orientation changes from an image to another but also because orientation trends change within images due to medium and high scales non stationarities especially in disordered materials such as C fiber and SL.

In figure 9, 𝑅𝐼𝑀𝑂𝐷(, ) is plotted as a function of distance for two angles ( = 0° and  = 90°) and for the four materials. The plots confirm the higher anisotropy of materials ReL and RL compared with SL and C Fiber . In all cases, orientation spatial differences are low for very short distances and increase to reach an asymptotic trend for long distances corresponding with the twist disorientation 𝛽 𝑀𝑂𝐷 . This asymptotic trend corresponds to high scale orientation variability. Differences can be observed in the rates of convergence towards the long range variability, the loss of coherence for ReL being faster than for RL. The twist disorientation 𝛽 𝑀𝑂𝐷 is also larger for ReL than for RL, respectively 10° and 8.5°, which is consistent with the orientation angles determined by selected area electron diffraction as shown in table 2. Discrepancies between angles =0° and =90° are also put into light, revealing that a preferential anisotropy is observed orthogonally to the fringes, especially for material RL. smaller than the values determined by XRD and the ratio 𝐿 𝑐 / 𝐿 𝑎 is larger for the MOD method than for XRD results. These points will be further discussed below.

Table 2: coherence lengths determined using the Mean Orientation Difference; average orientation angle 𝛽 𝑀𝑂𝐷 and coherence lengths determined by XRD for both PyCs. Provided uncertainties for the RIMOD method are the standard deviations estimated over 10 images.

While the curves in figure 9 differentiate between longitudinal and transverse coherence, it is also interesting to figure out how anisotropy varies according to angle . This can be done by plotting 𝑅𝐼𝑀𝑂𝐷 iso-levels using a radar representation. The graphs in figure 10 represent the distance ρ (in nanometers) at which feature 𝑅𝐼𝑀𝑂𝐷 reaches a value of θ degrees. θ is taken in the range 1° -9° for materials ReL and RL and in the range 0° -35° for materials SL and C Fiber . In the case of ReL and RL, a stronger spatial coherence is observed when  = 90° i.e. orthogonally to the fringes, which suggest that the stacks of near parallel fringes are non isotropic, their height being larger than their width. These iso-level curves also make easier the comparison between samples. For instance, it appears that ReL and RL have similar anisotropy at very short distances whatever the angle , but differ at larger distance, as RL seems to have bigger fringe stacks than ReL. This procedure allows obtaining a distribution of fringe sizes, which are shown at figure 12 for ReL and RL. It appears that the coherence length 𝐿 𝑎𝑀𝑂𝐷 is close to the mode of the distribution of fringe lengths while by XRD 𝐿 𝑎 values are shifted towards the largest fringe lengths. Two reasons could explain the larger values obtained by XRD. First, the CarbonXS program assumes a Gaussian distribution of stacking sizes whereas real domain size distributions are usually log-normal, leading to an unrealistic mean value. Second, the coherence lengths determined by full profile matching XRD methods [START_REF] Ruland | X-ray scattering of non-graphitic carbon: An improved method of evaluation[END_REF][START_REF] Faber | Investigation of the microstructure of disordered, non--graphitic carbons by an advanced analysis method for wide-angle Xray scattering[END_REF] are length weighted values (second moment the distribution divided by the first moment) [START_REF] Faber | Investigation of the microstructure of disordered, non--graphitic carbons by an advanced analysis method for wide-angle Xray scattering[END_REF], which are usually greater than the arithmetic mean; it can be seen in figure 2 that the length weighted average for the ReL PyC (𝐿 𝐿𝑊 -𝑅𝑒𝐿) is similar to L a (XRD) while the one for the RL PyC is slightly smaller than 𝐿 𝑎 (3 nm from the fringe length distribution and 4.1 nm from XRD).

Furthermore it must be underlined that XRD results are averaged over a large volume of material while TEM results were obtained on a much small amount and are thus less statistically significant. Thus the structural parameters 𝐿 𝑎𝑀𝑂𝐷 and 𝐿 𝑐𝑀𝑂𝐷 give an information averaged over the perfectly ordered and the defective areas (which are however weighted by the confidence factor) and are consequently smaller than the average fringe lengths. Besides, the large difference between the 𝐿 𝑐𝑀𝑂𝐷 /𝐿 𝑐𝑀𝑂𝐷 (~2) ratio and the 𝐿 𝑐(𝑋𝑅𝐷) /𝐿 𝑎(𝑋𝑅𝐷) ratio (~1) could be due to the projection effect. Indeed, the thin slice thickness is estimated to be in the range 10-30 nm and thus a large number of domains are superimposed along the thickness. Larger models (>10 nm) with known domain sizes will be needed to confirm this hypothesis.

Concerning materials C Fiber and SL, the rates of convergence are especially noteworthy. From figure 9 it can be observed that at short distances, the loss of coherence is slightly faster for SL (below  = 1 nm) but slows down afterwards and becomes slower than for material C Fiber that reaches its largest value (𝑅𝐼𝑀𝑂𝐷(, ) = 45.6°) for  = 7.5 𝑛𝑚, at a continuous rate; then, more surprisingly, for C Fiber , 𝑅𝐼𝑀𝑂𝐷 slightly decreases for larger distances ( 𝑅𝐼𝑀𝑂𝐷(, ) = 43.6° for  = 15 𝑛𝑚 ). This behavior can be explained by the specific microstructure of ex-PAN carbon fibers; they are made of a deep entanglement of folded sheets of layers, crumpled along the fiber axis, a microstructure fairly similar to the stacking of concentric cylinders artificially generated fig. 5. The curvature of the sheets are rather homogeneous as can be seen on the HRTEM image, fig. 6, which results in basic structural units with a mean curvature radius of roughly 8 nm, clearly visible in the local orientation map (fig. 7). The slight decrease of 𝑅𝐼𝑀𝑂𝐷 at larger distances may be due to the parallelism between neighboring BSUs which have a ribbon-like shape; however the oscillating feature shown in figure 5 is not observed in practice probably because the domains are not perfectly cylindrical and have various curvatures and sizes.

Conclusion

A new method of analysis of the HRTEM images of turbostratic carbons such as pyrocarbons or carbon fibers has been proposed. Contrary to standard procedures this method does not detect individual fringes but finds their local orientation and look at the spatial correlation between orientations. An analysis based on second-order statistics of the orientation field allows achieving a quantitative description of the nanostructure: coherence lengths 𝐿 𝑎𝑀𝑂𝐷 and 𝐿 𝑐𝑀𝑂𝐷 give the mean dimension of the BSU parallel and perpendicular to the fringes and 𝛽 𝑀𝑂𝐷 is a determination of the mean disorientation of fringes. A more thorough analysis of the mean orientation difference plots can give information such as the curvature of the USB and the correlation of orientations at large distance.

The method has been applied to the characterization of two high textured pyrocarbons. The results have been compared with some experimental data obtained by X-ray diffraction and electron diffraction. All techniques give consistent results and show that the rough laminar PyC exhibits larger coherent domains and is more textured at every length scale than the regenerative laminar PyC. Moreover the orientation field analysis provides a clear visualization of the defects, evidenced in the HRTEM images of high textured PyCs as line defects, perpendicular to the fringes and located at the domain boundaries. This method has already been used to validate some atomistic models of pyrocarbons [START_REF] Farbos | Nanoscale structure and texture of highly anisotropic pyrocarbons revisited with transmission electron microscopy, image processing, neutron diffraction and atomistic modelling[END_REF] and thus can help in the generation of more realistic models.

The application of the method on two less textured carbons (smooth laminar PyC and ex-PAN carbon fiber) showed that it can reveal information at different length scales: at short lengths scale (< 2 nm) the extent of domain sizes are revealed, while at larger length scale (> 4 nm) information such as the heterogeneity of the sample (SL PyC) or the presence of regularly spaced tubular structures (C Fiber ) are enlightened.

Future works will focus on two points: first, on a better recognition of the amorphous areas, holes and images of lattices planes perpendicular to the electron beam in order to extend the method to more complex images such as HRTEM images of soot or activated carbons and second, on the determination of domain sizes ( 𝐿 𝑎𝑀𝑂𝐷 , 𝐿 𝑐𝑀𝑂𝐷 ) on these less organized carbons.

  2 ) ∈ [0, 𝜋[× [0, 𝜋[ by using the function ∆(𝜃 1 , 𝜃 2 ) = min (|𝜃 1 -𝜃 2 |, 𝜋 -|𝜃 1 -𝜃 2 |). Based on ∆, various secondorder features were proposed to describe orientation dependencies in texture, among which the 𝑀𝑂𝐷 feature defined by:𝑀𝑂𝐷(𝐰) = 〈𝜂 𝐮 𝜂 𝐮+𝐰 ∆(𝜃 𝐮 , 𝜃 𝐮+𝐰 )〉 Ωwhere 〈•〉 Ω denotes spatial averaging over the whole image domain Ω. 𝑀𝑂𝐷 can be considered as a weighted measure of similarity between two orientations separated by spatial lag 𝐰 ∈ ℝ.We propose here to use a new rotation invariant formulation. Let (𝜌, ) be the polar representation of lag 𝐰. The rotation invariant formulation consists in considering the polar shift (𝜌, ) not in the image grid coordinate system but in a local system defined by the local orientation: 𝑅𝐼𝑀𝑂𝐷(𝜌, ) = 〈𝜂 𝐮 𝜂 𝐮+𝛅 ∆(𝜃 𝐮 , 𝜃 𝐮+𝛅 )〉 Ω with 𝛅 = (𝜌 cos(𝜃 𝐮 + ) , 𝜌 sin(𝜃 𝐮 + )).

Figure 1 .

 1 Figure 1. Polar representation of the spatial lag (𝜌, ) used in the case of the rotation invariant formulation 𝑅𝐼𝑀𝑂𝐷.

  (ii) a uniform [0,180°] random variable. Examples of such images are shown in figure 2. The full images have a size of 2048×2048 pixels and contain 2552 domains; the resolution was fixed to 0.037 nm/pixel i.e. the resolution of the experimental digitized HRTEM images used in this work. Thus the average grain area is 𝑆 𝑔 = 2.25 𝑛𝑚² and the average grain dimension 𝐷 = √ 𝑆 𝑞 = 1.5 𝑛𝑚. The two images shown in figure 2 only differ by the orientations of their fringes. The image in fig. 2a corresponds to a highly textured case while fig. 2d corresponds to an isotropic case. Images have been filtered with the radial frequency filter with 𝜌 0 = 1 𝑑 002 ≈ 2.9 nm -1 in order to smooth grain boundaries. Both orientation maps (Fig 2b and 2e) show that individual grain orientations are clearly identified while orientation confidence maps (Fig 2c and 2f) highlight the regions of disorder mainly located at grain boundaries.

Figure 2 :

 2 Figure 2: zoom of two artificial lattice fringe images (a,d) and corresponding local orientation (b,e) and orientation confidences (c,f). The parameters used for the simulation are: images of 2048×2048 pixels, 2552 grains, standard deviation of grain orientations  = 9° for picture (a) and uniform orientation distribution for picture (d). Both images were preprocessed with a band-pass radial filter. Resolution is 0.037 nm/pixel, 𝑑 002 = 0.34 𝑛𝑚.

Figure 3 :

 3 Figure 3: Plots of 𝑅𝐼𝑀𝑂𝐷(𝜌, 0°) and 𝑅𝐼𝑀𝑂𝐷(𝜌, 90°) for the two images in figure 2, high textured (left) and isotropic sample (right), dotted lines correspond to the parameters  MOD , L cMOD and L aMOD .

Figure 4 shows

 4 the corresponding 𝑅𝐼𝑀𝑂𝐷 diagram for columnar grains, it appears that curves at 0° and 90° are clearly different, leading to a ratio 𝐿 𝑐𝑀𝑂𝐷 /𝐿 𝑎𝑀𝑂𝐷 of 1.8 close to the shape anisotropy of 2 used to synthesize the image. Therefore 𝐿 𝑐𝑀𝑂𝐷 and 𝐿 𝑎𝑀𝑂𝐷 give information about domain sizes perpendicular and parallel to the layers.

Figure 4 :

 4 Figure 4: Mean orientation difference diagram corresponding to an image with columnar domains and zoom of the corresponding image and of the local orientation map. The size of the full image used for orientation statistics computing is 2048×2048 pixels (5742 nm²). It

Figure 5 :

 5 Figure 5: zoom of an artificial image made of a stacking of concentric circles, local orientation and Mean Orientation Difference diagram.

Figure 6 :

 6 Figure 6: raw HRTEM images of various carbons: upper left ReL, upper right RL, bottom left SL and bottom right C Fiber . Image size is 37.9nm×37.9nm. 9.5nm width zooms are shown as inserts.

Figure 7 :

 7 Figure 7: Local orientation for the four carbons (ReL, RL, SL and C Fiber , from top left to bottom right), orientation palette and zooms of the filtered HRTEM images.

Figure 8 :

 8 Figure 8: orientation coherences for the four carbons (ReL, RL, SL and C Fiber , from top left to bottom right).

Figure 9 :

 9 Figure 9: Plot of 𝑅𝐼𝑀𝑂𝐷(, ) as a function of distance , for angles =0° and 90°. Left: ReL and RL PyCs. Vertical dotted lines correspond to the coherence lengths 𝐿 𝑐𝑀𝑂𝐷 and 𝐿 𝑎𝑀𝑂𝐷 and horizontal dotted line to the average disorientation angle 𝛽 𝑀𝑂𝐷 . Right: SL PyC and C Fiber , the black dotted line corresponds to the maximum difference of orientation distance for the C Fiber .

Figure 10 :

 10 Figure 10: 𝑅𝐼𝑀𝑂𝐷 iso-level polar graphs computed on the four carbon images (ReL, RL, SL and C Fiber , from top left to bottom right). The graphs represent the distance ρ (in nanometers) at which feature 𝑅𝐼𝑀𝑂𝐷 reaches a value of θ degrees. θ is taken in the range 1° -9° for materials ReL and RL and in the range 0° -35° for materials SL and C Fiber .

Figure 11 :

 11 Figure 11: zoom of an HRTEM image for the ReL PyC (left) and detected fringes (right).

Figure 12 :

 12 Figure 12: Fringe length distribution determined for ReL (continuous line) and RL (dotted line). Abscissae corresponding to 𝐿 𝑎𝑀𝑂𝐷 , 𝐿 𝑎(XRD) are reported on the diagram as well as the length-weighted averages (L LW ), calculated from the fringe length distributions.

Figure 13 :

 13 Figure 13: Representation of a basic structural unit of a turbostratic carbon, according to [52], L 1 , L 2 and  (twist disorientation angle) obtained by high resolution transmission electron microscopy. 𝐿 𝑐 (average stack height) and 𝐿 𝑎 (average graphene layer extent) are determined by X-ray diffraction and are length-weighted averages.

Figure 14 :

 14 Figure 14: From left to right: slice of a model of pyrocarbon according to Leyssale et al. [8], projected potential and corresponding simulated HRTEM image.

Table 1 :

 1 determination of 𝐿 𝑐𝑀𝑂𝐷 , 𝐿 𝑎𝑀𝑂𝐷 and 𝛽 𝑀𝑂𝐷 from various artificial lattice fringe images.

	Image Size (pixels)	N	D (nm)	 °	𝐿 𝑐𝑀𝑂𝐷 (nm)	𝐿 𝑎𝑀𝑂𝐷 (nm)	𝐿 𝑐𝑀𝑂𝐷 𝐿 𝑎𝑀𝑂𝐷	𝛽 𝑀𝑂𝐷 °
	1024×1024 5742	0.50	9.0	0.54	0.47	1.14	13.2
	1024×1024 1346	1.03	9.0	1.06	0.97	1.09	11.7
	2048×2048 2552	1.50	9.0	1.48	1.44	1.03	11.2
	2048×2048 1436	2.00	9.0	2.06	2.05	1.01	11.0
	2048×2048	638	3.00	9.0	2.92	2.85	1.03	10.9
	2048×2048 2552	1.50	isotropic	1.55	1.54	1.01	44.9

  Coherence lengths 𝐿 𝑐𝑀𝑂𝐷 and 𝐿 𝑎𝑀𝑂𝐷 are determined for both PyCs, 𝐿 𝑐𝑀𝑂𝐷 = 0.99 𝑛𝑚 and 𝐿 𝑎𝑀𝑂𝐷 = 0.51 𝑛𝑚 for ReL and 𝐿 𝑐𝑀𝑂𝐷 = 1.27 𝑛𝑚 and 𝐿 𝑎𝑀𝑂𝐷 = 0.55 𝑛𝑚 for RL. These results show that domains are larger for the RL PyC with a more columnar shape as compared with the ReL PyC. It is of interest to compare these results with the coherence lengths 𝐿 𝑐 and 𝐿 𝑎 determined by X-ray diffraction; results are given table 2 (Standard deviations, given in table2have been estimated through the analysis of a set of ten HRTEM images, thus taking into account the spatial variation inside the sample). It appears that coherence lengths obtained by XRD and by MOD are larger for RL than for ReL, however 𝐿 𝑐 and 𝐿 𝑎 obtained by MOD are

Acknowledgements

Funding from the French Agence Nationale de la Recherche through the PyroMaN project (contract ANR-2010-BLAN-929) and from the Carnot Materials and systems Institute of Bordeaux (MIB) are gratefully acknowledged. The authors are grateful to Professor Jeff Dahn for kindly providing the CarbonXS program.


Javier Rojas Balderrama

Matthieu Simonin

Cédric Tedeschi

A Chemistry-inspired Programming Model for Adaptive Decentralised Workflows

Keywords: Scientific workflows, adaptiveness, Chemical Programming Model Workflows scientifiques, adapation, modèle de programmation chimique

In this paper, we devise a chemistry-inspired programming model for the decentralised execution of scientific workflows, with the possibility of dynamically adapting its shape when its initial specification fails to reach the user's requirements or simply to run due to external conditions. We describe a decentralised architecture to support the model and cover its implementation in the GinFlow software prototype.

Introduction

In contrast to business workflows, whose outcome is quite predictable, the output of scientific workflows (and their potential intermediate data) is rather hard to characterise prior to execution. This uncertainty comes together with the scientist's requirement to be able to explore different scenarios during computation time, due to some (possibly unsatisfactory) intermediary result. Therefore, scientific workflows are meant to support fields where variability is still largely present during experimentation.

Consequently, a workflow specification language should be able to deal with alternative workflow scenarios which will effectively get enabled only if some specific event is detected, such as a failure or the non-satisfaction of some specific property of some partial data results. This would provide a new chance to obtain meaningful results without having to start over the whole workflow enactment. At run time, this alternative workflow should be automatically triggered by the workflow engine, and in a transparent way for the user.

Another common need is scalability support. Scientists cannot afford being limited by workflow engines that may not be able to cope with workflows requiring a high overhead due to coordination at run time. Most workflow manager systems are centralised, and thus suffer from all the limitations inherent to it. However, building a decentralised engine raises several challenges related to the possibility to coordinate tasks together in a transparent way for the user.

In this paper, we tackle the general problems of finding 1) adequate programming abstractions and 2) system architectures for the development of a decentralised workflow management system with adaptiveness capabilities. Concretely, we devise the design, programming model and implementation of a decentralised workflow engine, which makes it possible to update the workflow while it is running, provided the alternative definitions to replace unsatisfactory/failed part of the workflow. Contributions. We formulate a set of programming abstractions using a declarative chemistryinspired model. It provides an elegant and concise way of specifying self-adaptive programs through the higher order (rules injecting new rules at run time) to define a scientific workflow, the data composition strategies, and the workflow variants declared in case of exceptions. We define two families of definitions, one for abstract workflows, to use at the user-level; another for concrete workflows, to instantiate the abstract workflow specification and enable the adaptation rules. Finally, we review the implementation and the software technologies it relies on. Outline. The rest of the paper is organised as follows. The programming abstractions for an adaptive decentralised workflow are detailed and exemplified in Section 2. The architecture and implementation is discussed in Section 3. The closing sections deal with related works, and conclusions.

Programming model

Artificial chemistries [START_REF] Dittrich | Artificial chemistries-A review[END_REF] provide a programming style based on the chemical metaphor. They have been shown relevant for the specification of autonomous systems. In such programming models, data are envisioned as molecules floating in a chemical solution, and reacting according to reaction rules (i.e., the program) to produce new molecules (the resulting data). Reactions are conditional, and take place between a (multi)set of molecules satisfying a reaction condition. This process continues until no reactions can be performed anymore so the solution is inert. Reactions take place in an implicitly parallel and autonomous way, and in a non-deterministic order. Formally, the solution is a multiset of molecules, and the possible reactions between molecules are a set of rewriting rules on this multiset. The multiset is the unique data structure, and acts similarly to an address space through which processors can cooperate. In this section, we define a set of programming abstractions for adaptive decentralised workflows. We present HOCL flow , an extension of the HOCL language. HOCL flow enhances the original specification with 1) the support for some data structuring, and 2) macros to improve the clarity and conciseness of the rules.

Higher-Order Chemical Language (HOCL).

HOCL is a higher-order, rule-based declarative language inspired by a chemical metaphor and proceeding by state rewriting [START_REF] Banâtre | Generalised multisets for chemical programming[END_REF]. In HOCL, a program is a solution (i.e., a multiset) of atoms A 1 , . . . , A n . An atom can be a constant (integers, booleans, etc); a tuple of atoms e 1 : e 2 : • • • : e n ; a subsolution B 1 , . . . , B m of atoms, or a reaction rule. A multiset of atoms is called a molecule. Note that the solution of a program is the greatest molecule of this program. A possible reaction rule between elements of the solution is written replace-one P by M if C, where P denotes the pattern to be satisfied by the molecules and C the condition they must satisfy. If that is the case, the reaction consumes the rule and the matched molecule, produces M , and inject it into the solution. A replace-one rule is one-shot: it disappears when it reacts. Its variant replace P by M if C is n-shot: it is not consumed when it reacts. For instance, consider a solution which calculates the maximum value out of a given set of numbers. The example below illustrates the expressiveness and higher order of HOCL, where reactions consume and/or produce other reaction rules. Let us start with the following simple program:

let max = replace x, y by x if x ≥ y in 2, 3, 5, 8, 9, max
The max rule consumes two integers x and y when x ≥ y and replaces them by x. Initially, several such reactions are possible in the provided multiset, max can use any couple of integers satisfying the condition: 2 and 3, 2 and 5, 8 and 9, etc. We may introduce a higher-order rule to only obtain the result in the final solution. Note that it is responsible to delete the max rule once the solution only contains the highest integer value. This introduces the need for sequentiality of events because we need to wait that all possible application of max to take place before deleting it. Within the chemical model, the sequentiality is achieved through subsolutions. A rule has to Inria wait for its inertia to access a subsolution. In our example, this leads to the encapsulation of the solution and the following program: 2, 3, 5, 8, 9, max , replace-one max = m, ω by ω

The m variable matches a rule named max, and ω matches all the remaining elements. One possible execution scenario within the subsolution is the following (2 and 8 on one hand, and 3 and 5 on the other hand react first, producing the intermediate state): 2, 3, 5, 8, 9, max → * 3, 5, 9, max → * 9, max Once the inertia is reached within the subsolution, the one-shot rule is triggered, extracting the result:

9, max , replace-one max = m, ω by ω → 9 This fine-grain example shows that HOCL provides the ability to express autonomic coordination of rules (without necessitating a centralised control). The current state of a computation is represented by the solution, which constitutes an information system by itself. The multiset becomes a shared space providing the information required for dynamic coordination, a suited requirement in a decentralised workflow engine.

HOCL flow

We here define HOCL flow , an extension of HOCL. The original specification supports natively very limited data structuring, through subsolutions -hierarchical structuring of unstructured subsets of elements-and tuples-ordered molecules. Tuples however can only be manipulated through pattern matching which sometimes leads to complex patterns in the left part of rules definition. We introduce lists and basic manipulation functions to improve data manipulation needed between the tasks of a workflow (see Section 2.3). Given a list and an element e, we define the following:

• () is a empty list • list(e) creates a list containing the element e

• cons(e, ) appends e to the list • first( ) returns the first element of • rest( ) returns a list whose value is without its first element • nth( ) returns the n-th element of We also define the macro construct inject. HOCL flow is designed to make rules easier to read or write. In other words, it provides syntactic sugar on top of HOCL. This macro adds content into a solution without having to repeat molecules that appear in both left and right parts of the rule:

inject M ≡ replace-one ω? by M, ω When a catalyst T is needed (i.e., an extra molecule participating in the rule activation without being consumed or modified), the with clause can be added to the inject macro along with an optional condition C:

inject M with T if C ≡ replace-one T, ω? by M, T, ω if C
where C requires the contents of T to activate the rule (e.g., the if clause on some element in T ). Nonetheless, the with-clause can be used without any conditional when parts of contents of T are simply present in M .

Tasks, bindings and dependencies

There are two fundamental information without which a workflow would not exist and which are, consequently, required by any workflow runtime system. Firstly, the binding to actual services performing the tasks; and secondly, the links connecting these tasks together (e.g., the successors and predecessors for the task). Let us take a simple workflow example and review its HOCL flow representation, as illustrated by Figure 1.

A workflow is composed of as many subsolutions as tasks in a DAG definition. Each subsolution contains binding information and dependencies for one task. The atom srv declares the binding to the actual service s i . The molecule in : contains the parameters required to invoke it, and molecules src : and dst : the source and destination dependencies, respectively. For instance, T 2 requires the results of T 1 as input, and its results will be sent to T 4 .

A workflow definition needs to get set into action, since the above HOCL code is a mere description of what to do (not how to do it). Any workflow engine needs to 1) call services, and 2) transfer data between them. Rules are re-writers, as previously explained, so they can change the workflow state expressed by the HOCL flow multiset. The application of such rules allow us to reflect the completion of a task or the availability of partial results. Hence, rules should be generic to process any HOCL flow workflow definition. In order to achieve those two requirements, three generic workflow rules (gw_ * ) are defined below:

let gw_setup = replace-one src : , in : ωin? (1)
by src : , par : par = list(ωin), res : () let gw_call = replace-one src : , srv : si, par : par, res : res

by src : , srv : si, par : rest( par), res : cons(invoke(si, first( par)), res) if par = () let gw_pass = replace Ti : par : par, res : res, dst : Tj, ωdst? , ωi? ,

Tj : src : Ti, ωsrc? , in : ωin? , ωj?

by Ti : par : par, res : res, dst : ωdst , ωi , Tj : src : ωsrc , in : res, ωin , ωj if par = ()

Rule gw_setup detects the readiness of the task (i.e., all dependencies have been satisfied) through the emptiness of the src : molecule. Upon application, the rule moves the content from the molecule in : to the atom par as a list of parameters. It also creates the atom res as an empty list which, later, will contain the results of the service invocation. These molecules altogether trigger the gw_call which invokes the service implementing the task. After service completion, the result is inserted in res. This rule is applied as long as elements are available in par. At first sight, in : and par seem to be redundant, however, in : is a container for results received from sources, whereas par is a list buffer of parameters for the invocation of services. Section 2.3 describes their different purpose and usage in greater detail. In practice, rules gw_setup and gw_call are used together. They are internal in the sense that they have to be present in each subsolution (see Figure 2).

Rule gw_pass is responsible for transferring results from a source to forthcoming destinations. It is global to the workflow because its scope is not limited to one subsolution. Its scope spans at least two services (one source and one destination). gw_pass is triggered after the result has been obtained and placed in res. It moves the resulting value from the source as inputs to each declared destination. Then it updates the sources and destinations to reflect that the dependency is satisfied. This rule may be triggered as many times as necessary according to all dependencies defined in molecules src : and dst : . These generic rules constitute a high-level programmatic way to enact a workflow expressed with HOCL flow . In means that an abstract workflow is changed into a concrete one by adding these rules. Once the concrete workflow is obtained, it can be run by the HOCL interpreter. The concrete workflow corresponding to the workflow in Figure 1 is provided in Figure 2.

Data composition strategies

It is common that a task does not receive a simple result from its sources but a data set. How to combine the elements of these sets into distinct inputs to the service realising the tasks relies on composition strategies. In this section, we introduce the HOCL flow mechanisms to process some of the most common strategies described in literature [START_REF] Oinn | Taverna: lessons in creating a workflow environment for the life sciences[END_REF], namely the dot and cross products, and the input filtering. In the following, let us assume tasks have results provided as lists, for example, res i = (a 0 , a 1 , . . . , a m-1 ), and res j = (b 0 , b 1 , . . . , b n-1 ).

Dot product.

The dot product is the set of ordered pairs (a k , b k ) for all k such that 0 ≤ k < min(m, n) from the resulting combination of the elements of res i and res j represented as lists i and j respectively. Each pair serves as an independent input for the invocation of the service supporting the task. In HOCL flow , the dot product (dp_ * ) can be defined through rules 4 and 5.

let dp_cons = replace src : dot , in : i, j , par : par (4)

by src : dot , in : rest( i), rest( j ) ,

par : cons((first( i) : first( j )), par) if i = () ∧ j = () let dp_setup = replace-one src : dot , in : i, j (5)
by src : , res : ()

if i = () ∨ j = ()
At run time, the rule dp_cons is applied as soon as the required molecules are available. It creates the pairs for the service invocation in the par atom. This allows the interpreter to apply the rule gw_call, as many times as necessary-once for each pair. Consequently, the rule dp_setup can be applied when processed molecules are empty, paving the way for the application of the gw_call rule with the created pairs, again once for each pair.

Cross Product.

The cross product is the set of all ordered pairs (a k , b l ) with 0 ≤ k < m and 0 ≤ l < n. Each pair serves as an independent input for an invocation of the service supporting the task. In HOCL flow , the cross product (cp_ * ) can be defined through rules 6, 7, 8, and 9.

let cp_init = inject in : backup = j with in : i, j (6)
let cp_cons = replace src : cross , in : i, j , par : par

by src : cross , in : i, rest( j ) , par : cons((first( i) : first( j )), par) if j = ()

let cp_next = replace src : cross , in : i, j , in : backup (8)
by src : cross , in : rest( i), backup , in : backup

if i = () ∧ j = () let cp_setup = replace-one src : cross , in : i, j , in : backup (9)
by src : , res : ()

if i = ()
The cross product is handled in a way similar to the dot product. Firstly, the cp_init rule creates a backup copy of the j list, by injecting a in : molecule. This copy is used after each element in the list i has been combined with all of the elements in j , to refill in : so it can be combined with the next element in i . Then, cp_cons iterates through j while the element of i is fixed, removing the element from j at each iteration. Once j is empty, the rule cp_next can be applied to refill j and remove the element of i that was just processed. The rule cp_setup has the same role as the homologous rule dp_setup.

Filtering.

A composite result may be composed of several parts having different purposes. In workflows, this is commonly dealt with through data filtering: when the result reaches one of its destination, it goes through a filter so only the relevant part in the source is selected. The filtering (df_ * ) can be defined through rules 10 and 11.

let df_init = replace-one src : filter : (Ti : j), ωsrc? , in : i, ωin? At the abstract workflow level, expressing a filter requires adding a molecule filter : (T i , j) on src : . The subsolution T i denotes the origin of the composite result, and j gives the rank of the element needed in i . Each result can be filtered at several positions, with as many filter : i , * ) molecules as necessary to filter i . The df_init rule isolates the i list provided a molecule filter : (T i , * ) is found. Note that even if several molecules filter : (T i , * ) are present, performing this isolation once is enough. Then, the rule df_step can be triggered as many times as necessary to actually keep only the requested parts.

Example.

Let us consider the simple workflow described in Figure 3 to illustrate the building blocks of our programming model presented above. This workflow has the same structure as the workflow depicted previously in Figure 1. Its difference stands in the sense that it locally pre-processes data prior to some of services' invocation. Specifically, it applies a filter on the result of T 1 before submitting it to Service s 2 . A dot product is also applied on the results of T 2 and T 3 before invoking T 4 . Note that only the relevant generic rules were inserted. The details of the generation of the concrete workflow are described in Section 3. We detail next the execution of the workflow focusing on the rewriting of the multiset. For the sake of clarity, the generic rules are omitted in the multiset, so only the abstract workflow is represented. Initially, only the gw_setup rule is enabled, in subsolution T 1 . It is triggered to prepare the service invocation. An atom par is created, rewriting the subsolution T 1 :

T 1 : src : , srv : s 1 , par : (input), dst : T 2 , T 3 Again, only the gw_call rule can be applied. Triggering it will actually call the service and inject the result in Subsolution T 1 :

T 1 : src : , srv : s 1 , par : (), res : (res 1 ), dst : T 2 , T 3 Heretofore, subsolutions other than T 1 have been inert: no combinations of molecules in them could match one rule's left part. The outside rule gw_pass is now enabled; the results of T 1 are available and they can be sent to its destinations, namely T 2 and T 3 . Note that the gw_pass rule is triggered twice, once for each destination of T 1 , adding to in : the respective res 1 . At this point, subsolution T 1 is inert definitively. All these subsolutions are then updated, reflecting the transfer and leading to the following multiset: T1 : src : , srv : s1, par : (), res : (res1), dst : T2 : src : filter : (T1 : 3) , srv : s2, in : res1 , dst : T4 T3 : src : , srv : s3, in : res1 , dst : T4 T4 : src : dot, T2, T3 , srv : s4, in : , dst :

RR n°8691

It is now possible to process T 2 and T 3 in parallel. In T 3 , rules gw_setup and gw_call can be triggered, one after another. However, due to the presence of filter : (T 1 : 3) inside src : , Rule gw_setup is not enabled in T 2 because it requires an empty molecule src : . That same atom will first trigger the df_init and df_step rules, extracting the filtered part of the result. At the same time, src : is flushed, and a par atom is created. Therefore, the subsolution T 2 is now ready for the application of the rules gw_setup and gw_call: Now dp_cons can be applied and the pairs are created by its repeated application. Let us assume, for instance, that the results of subsolutions T 2 and T 3 are made of two elements each, denoted res 21 , res 22 and res 31 , res 32 respectively. These results in the subsolution T 4 are defined as following:

T 4 : src : dot , srv : s 4 , in : , par : ((res 21 : res 31 ), (res 22 : res 32 )), dst :

Once dp_setup has been applied, cleaning src : , gw_call is triggered for each pair of par, invoking s 4 twice, and obtaining the two final results. This leads to the final state of all workflow and an inert multiset.

Adaptiveness

Exception-handling support during workflow execution is typically managed at infrastructure level. Autonomic systems adapt to exceptional situations while hiding their intrinsic complexity to users. In scientific workflows, those exceptions (e.g., unavailability, software error or middleware failures) can be detected by monitoring the execution status of the services and providing corrective actions using pre-defined policies. An alternative approach in programmatic environments may be to declare the adaptiveness in an explicit manner. Such a programmatic adaptiveness provides mechanisms to adapt the workflow itself when facing changing conditions while letting users the possibility to express the way the workflow enactment evolves. Similarly to the basic workflow treatment, the user-defined abstract adaptive workflow needs to be extended becoming a concrete workflow as executable HOCL code.

We identify four steps to implement workflow adaptiveness as an extension of the set of rules defining the workflow engine. These are: 1) the declaration of alternatives to replace failing services 2) the identification of the exception, 3) the propagation of the exception to the concerned dependent services, and 4) replacement of concerned services with their equivalent, replacement services. The user is responsible for the inclusion of an alternative scenario in the workflow definition. An alternative workflow excerpt must be either a simple service or a composite service, but with only one exit route. This condition is set to avoid incoherent states during the enactment. A workflow example illustrating the adaptiveness is given in Figure 4.

The adaptation strategy is declared in an abstract workflow by defining the subsolution(s) associated with the alternative service(s). The main idea is to rely on an inline rule that, combined with the replacing service(s), starts the adaptive process. This rule includes the atoms trigger and bond defining the replacement, the new bindings between services, and the exception condition. It is applied as soon as the propagation of an exception is detected. Atoms trigger and bond can be seen as special meta-clauses that will be replaced, during the generation of the concrete workflow out of the abstract one, with parametric rules to be enacted if the alternate workflow needs to get started. The trigger atom refers to the alternative services that may be used in case of another service raises an exception. The bond atom should specify where the adaptiveness needs to get propagated, so as to link the new dependencies (i.e., recreating the chemistry bonds between subsolutions). The inline rule has to declare as many atoms of type bond as dependencies are declared in the sources and destinations of the related subsolutions. This is all the information users need to define.

Let us review the workflow adaptiveness generic or generated rules (wa_ * ) that will be added to the user-defined workflow so as to enable the exception handling with high-level definition prior to execution. They are defined below:

let wa_init = inject dst : (12)
let wa_pass = replace Ti : par : par, res : res, dst : Tj, ωdst? , dst : ω dst ? , ωi? , Tj : src : Ti, ωsrc? , in : ωin? , ωj?

by Ti : par : par, res : res, dst : ωdst , dst : Tj, ω dst , ωi , Tj : src : ωsrc , in : res, ωin , ωj

if par = () ∧ failed / ∈ res (13)
let wa_update-src ij = replace-one src : Ti, ωsrc? by src : Tj, ωsrc

let wa_update-dst ij = replace-one dst : ωdst? , dst : Ti, ω dst ?

by dst : Tj, ωdst , dst :

ω dst let wa_trigger-dst kij = replace-one T k : ω k ? (16)
by

T k : wa_update-src ij , ω k let wa_trigger-src kij = replace-one T k : ω k ? (17)
by

T k : wa_update-dst ij , ω k
The wa_init rule injects, into each subsolution, the dst : molecule. This molecule acts as a supplementary holder of sources identified as having an alternative service in case of exception. Then the wa_pass rule, an extended version of its gw_pass analogous, not only transfers results from sources to destinations when the service invocation is achieved properly, but also updates all dependencies associated to the service. It replaces the service with the alternative declared in the excerpt associated to the workflow adaptiveness.

Rules wa_update- * and wa_trigger- * represent the core instrumentation of the workflow adaptiveness. Unlike all previous rules, this set of rules is parametric; they incorporate, in their definition, variables indexed i, j, k to refer to the involved subsolutions. They are generated specifically for each workflow based on the definition of the abstract workflow. Rules wa_update- * update the dependency bindings of the subsolutions changing the contents of sources and destinations. The wa_update-src ij rule updates the source T i with the alternative T j in all subsolutions containing this dependency. In the same way, wa_update-dst ij replaces the old references of the destination T i with the value T j .

Finally, rules wa_trigger- * , which take advantage of the high-order properties, adapt the workflow definition, enabling the new subsolution as a proper part of the workflow. The wa_trigger-dst kij rule amends the contents of subsolution T i containing T k as part of their destinations with the new reference T j . This process resets the inert status of already invoked services, by forcing to transfer its results to their newly added destinations, as defined in Rule wa_pass. In a similar way, the wa_trigger-src kij rule removes the references to the services which raised an exception, replacing them by the alternative subsolution. Following the example of Figure 4, the resulting concrete workflow will include the following excerpt: T4 : src : T1 , srv : s4, in : , dst : T3 , wa_init, gw_setup, gw_call inject wa_trigger-dst 324 , wa_trigger-src 124 with T2 : res : (ωres), ω2 if ωres = failed At time of concretisation, the subsolution and the inline adaptive rule will be used to generate the following set of rules: let wa_trigger-dst 324 = replace-one T3 : ω3? by T3 : wa_update-src 24 , ω3 let wa_trigger-src 124 = replace-one T1 : ω1? by T1 : wa_update-dst 24 , ω1 let wa_update-src 24 = replace-one src : T2, ωsrc? by src : T4, ωsrc let wa_update-dst 24 = replace-one dst : ωdst? , dst : T2, ω dst ? by dst : T4, ωdst , dst : ω dst In general, a failing service may be replaced with another equivalent service (or set of services). This replacement does not require any extra runtime procedures at the system level. The parametric rules define all the actions to be taken (provided an HOCL flow interpreter). The final definition of the adaptive workflow does not modify the original workflow definition. All the complexity of the enactment instrumentation is also hidden from the user.

Architecture and implementation

In this section, we describe the design principles of GinFlow, the software prototype developed to implement the proposed programming model in a decentralised fashion. As we explained previously, decentralisation is one of our main motivations. The architecture proposed takes its roots in [START_REF] Fernández | Rule-driven service coordination middleware for scientific applications[END_REF], and relies on a shared-space based coordination model, as depicted in Figure 5. We assume each service is deployed and encapsulated into a Service Agent (SA). SAs communicate Inria 6: Ginflow decentralised architecture through the shared space, which contains a description of (the state of) the workflow. This description is assumed to reflect the progress of the workflow, in particular completed tasks, and available input data. Each SA can rely on the space to act. The service agents execute one task, and take part in the coordination by reading and writing the shared space. Doing so, they notify when agents of subsequent executions can start (provided other incoming dependencies have also been satisfied).

The GinFlow prototype

Figure 6 shows a logical view of the GinFlow architecture. As mentioned before, the initial workflow is described as an HOCL program. The shared space is implemented through the multiset containing it. To implement the SAs, GinFlow depends on a set of HOCL interpretors acting as co-engines and a set of workers responsible for the actual execution of services. GinFlow may deploy as many co-engines as tasks in the workflow. Each co-engine concurrently read the multiset to collect (and copy locally) the status of the workflow. Each embedded HOCL engine can therefore detect if some process is needed.

For instance, when a dot product is computed on a co-engine, it iterates over the inputs as shown in Section 2.3 and produces a list of parameters. Then, it triggers the execution of the service, wrapped in one or several workers in parallel, depending on the size of the input parameters list. In other words, each time the gw_call rule is triggered, a new worker receives work. Finally, the co-engine collects the results and pushes them back, along with the adapted HOCL description of the new workflow status, into the shared multiset. The co-engines lay at the interface between the HOCL domain and the services. They act as local coordinators for the execution of the HOCL service described in the multiset through the workers. These workers are generic in the sense of execution because they can call different implementations (SOAP, REST, binaries, etc). They are handled through a common abstraction layer. Moreover they are lightweight, stateless and fully decoupled from the other components of the architecture. In the current implementation, the HOCL core language, and the messaging protocol inherent to the distributed nature of the implementation are independent modules. GinFlow relies on the ActiveMQ middleware to implement the communications between the multiset and the HOCL engines (i.e., the gw_pass rule is implemented through ActiveMQ queues).

Concrete workflow generation and enactment

Section 2.1 describes how abstract workflows are translated into concrete workflows prior to execution. The user only has to specify the source and destination for each service. The preprocessing phase will generate what is needed for the workflow to run correctly (passing molecules from one service to another, iterating over the inputs, and so on). Thus the compilation of an HOCL flow workflow follows a pipeline of operations that consists in 1) parsing the abstract description and adding the generic rules allowing service invocation, data pre-processing and transfer, 2) generating the concrete rules from meta-molecules (e.g bond and trigger), and 3) issuing the resulting workflow to the HOCL flow compiler.

The pre-processing steps (1) and ( 2) parse the input abstract workflow and add the relevant rules to the workflow. Following the same principle, specific rules supporting the adaptiveness are hidden from the user (see Figure 4). In order to implement this process, we rely on a mechanism similar to the one which enables the molecule transfer between services (using the gw_pass rule). Such communications between HOCL entities make use of the publish/subscribe capabilities of the ActiveMQ broker. Figure 7 illustrates the enactement of the adaptiveness process. When a service/co-engine detects a failure, it injects the updates rules in the subsolution of the common output T out and in the sources T in at runtime. As discussed in Section 2.4, the workflow will resume automatically.

Related work

In contrast to classical workflow manager systems, approaches using the chemical metaphor [START_REF] Caeiro | A chemical model for dynamic workflow coordination[END_REF][START_REF] Di Napoli | Adaptive instantiation of service workflows using a chemical approach[END_REF][START_REF] Németh | Distributed workflow coordination: molecules and reactions[END_REF] envisage the workflow execution as an autonomous process evolving in time according to the requirements and dependencies without bounding to any preset constraint. For instance, the work presented in [START_REF] Németh | Workflow enactment based on a chemical methaphor[END_REF] represents the enactment as an abstract coordination modelling which identify different levels of the workflow definition process: an abstract workflow that expresses the logic, and a concrete workflow where logical entities are assigned to resources for execution. They propose a uniform declarative formalism to describe workflows and coordination strategies. However, they do not exploit higher-order constructs nor envisage full dynamicity at runtime to support complex mechanisms such as fault-tolerance.

Another chemical model [START_REF] Caeiro | A chemical model for dynamic workflow coordination[END_REF] defines a general workflow notation, including common workflow constructs, separating data and control flow. This formalisation pays special attention to the resource instrumentation through its identifiers to ensure dynamic coordination. However, they do not address a decentralised execution. In the same direction, a service-based application Inria Our architectural approach takes its roots in the work presented in [START_REF] Fernández | Rule-driven service coordination middleware for scientific applications[END_REF]. Nevertheless, that previous work does not address data iteration strategies, and neither does adaptiveness. The GinFlow prototype also leverages workers for task parallelism, thus boosting its scalability in contrast to the implementation presented in [START_REF] Fernández | Rule-driven service coordination middleware for scientific applications[END_REF].

Most of workflows manager systems initiatives achieve the required enactment flexibility by means of the infrastructure (e.g., re-submissions strategies [START_REF] Lingrand | Modeling user submission strategies on production grids[END_REF], pilot jobs [START_REF] Casajus | DIRAC pilot framework and the DIRAC workload management system[END_REF]). On the other hand, Tolosana-Calasanz et al. [START_REF] Tolosana-Calasanz | Adaptive exception handling for scientific workflows[END_REF] propose an adaptive exception handling for scientific workflows at definition level. Their work proposes two patterns to manage the exception handling based on the Reference Nets-within-Nets formalism: propagation and replacement. In spite of mechanisms for dynamically adapting the workflow structure at runtime without having to be aware of the underlying infrastructure, the resulting representation with their reference model suggests a quite complex workflow definition, where the original scenario and the alternative path are mixed.

Conclusions

In this paper, we have presented a set of programming abstractions based on a chemical metaphor to express data-driven workflows. We introduced an extension of the HOCL language aiming at easing the development of such programs. We have shown how to, out of the mere description of the workflow, it can get enabled in a decentralised environment. Also, we devised a set of abstractions (and their implementation) to automate exception-handling, guided by the user at design time. The GinFlow prototype is currently under development and testing. In parallel to this work, the GinFlow prototype is currently integrated to the TIGRES workflow tools [START_REF] Balderrama | Combining Workflow Templates with a Shared Space-based Execution Model[END_REF].

T 1 :T 2 :T 3 :T 4 :

 1234 src : , srv : s 1 , in : input , dst : T 2 , T 3 , src : T 1 , srv : s 2 , in : , dst : T 4 , src : T 1 , srv : s 3 , in : , dst : T 4 , src : T 2 , T 3 , srv : s 4 , in : , dst :

Figure 1 :

 1 Figure 1: Workflow DAG (left), HOCL flow definition (right)

Inria

T 1 : 03 T 2 : 04 T 3 : 05 T 4 :

 1032043054 src : , srv : s 1 , in : input , dst : T 2 , T 3 , gw_setup, gw_call , 4.src : T 1 , srv : s 2 , in : , dst : T 4 , gw_setup, gw_call , 4.src : T 1 , srv : s 3 , in : , dst : T 4 , gw_setup, gw_call , 4.src : T 2 , T 3 , srv : s 4 , in : , dst : , gw_setup, gw_call

Figure 2 :

 2 Figure 2: Concrete workflow (abstract workflow + generic workflow rules)

Figure 3 :

 3 Figure 3: Concrete workflow with data pre-processing

T 2 :

 2 src : , srv : s 2 , in : res 13 , dst : T 4Once gw_setup and gw_call have been applied within subsolutions T 2 and T 3 , gw_pass can transfer results into in :in T 4 , as previously from T 1 to T 2 and T 3 . Subsolutions T 2 and T 3 are now definitively inert: T1 : src : , srv : s1, par : (), res : (res1), dst : T2 : src : , srv : s2, par : (), res : (res2), dst : T3 : src : , srv : s3, par : (), res : (res3), dst : T4 : src : dot , srv : s4, in : (res21 : res22), (res31 : res32), dst :

Inria

T 1 :T 2 :T 3 : 04 T 4 :

 123044 src : , srv : s 1 , in : input , dst : T 2 , src : T 1 , srv : s 2 , in : , dst : T 3 , src : T 2 , srv : s 3 , in : , dst : , 7.src : T 1 , srv : s 4 , in : , dst : T 3 , 7.05 inject trigger : T 4 , bond : (T 1 : T 3 ) 7.06 with T 2 : res : (ω res ), ω 2 if ω res = failed

Figure 4 :

 4 Figure 4: Adaptiveness: Schema (left) and abstract workflow (right)

Figure 5 :

 5 Figure 5: A decentralised workflow execution environment

Figure

  Figure 6: Ginflow decentralised architecture

Figure 7 :

 7 Figure 7: Adaptiveness implementation

RR n°8691

Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399