Gaüzère Benoit

Brun Luc

Villemin Didier

Two New Graphs Kernels in Chemoinformatics

Chemoinformatics is a well established research field concerned with the discovery of molecule's properties through informational techniques. Computer science's research fields mainly concerned by chemoinformatics are machine learning and graph theory. From this point of view, graph kernels provide a nice framework combining machine learning and graph theory techniques. Such kernels prove their efficiency on several chemoinformatics problems and this paper presents two new graph kernels applied to regression and classification problems. The first kernel is based on the notion of edit distance while the second is based on subtrees enumeration. The design of this last kernel is based on a variable selection step in order to obtain kernels defined on parsimonious sets of patterns. Performances of both kernels are investigated through experiments.

Introduction

Chemoinformatics aims to predict or analyse molecule's properties through informational techniques. One of the major principle in this research field is the similarity principle, which states that two structurally similar molecules should have similar activities and properties. The structure of a molecule is naturally encoded by a labeled graph G = (V, E, µ, ν), where the unlabeled graph (V, E) encodes the structure of the molecule while µ maps each vertex to an atom's label and ν characterizes a type of bond between two atoms (single, double, triple or aromatic).

1

A first family of methods introduced within the Quantitative Structure-Activity Relationship (QSAR) field is based on the correlation between molecule's descriptors and molecule's properties (e.g. molecule's boiling point). Vectors of descriptors may be defined from structural information (Cherqaoui and Villemin, 1994), physical properties or biological activities and may be used within any statistical machine learning algorithm to predict molecule's properties. Such a scheme allows to benefit from the large set of tools available within the statistical machine learning framework. However, the definition of a vector from a molecule, i.e. a graph, induces a loss of information. Moreover, for each application, the definition of a vectorial description of each molecule remains heuristic.

A slightly different approach is based on graph embedding. Within this framework, a vectorial description of a graph is automatically build from its encoding. A vectorial embedding of a graph may, for example, be defined by associating to each graph its distance to a set of graph prototypes [START_REF] Emms | Graph embedding using quantum commute times[END_REF]. Another graph embedding approach is based on the spectral analysis of graphs [START_REF] Caelli | An eigenspace projection clustering method for inexact graph matching[END_REF][START_REF] Luo | A spectral approach to learning structural variations in graphs[END_REF]. In this last case, the embedding is deduced from the analysis of the eigenvectors and eigenvalues of the adjacency matrix. This rich framework is closely related to several graph embedding methods such as the ones based on random and quantic walks [START_REF] Berchenko | Graph embedding through random walk for shortest paths problems[END_REF][START_REF] Emms | Graph embedding using quantum commute times[END_REF] or the ones based on heat kernel equations (Elghawalby andHancock, 2008, 2009).

A third family of methods, based on graph theory, may be decomposed in two subfamilies. The first subfamily [START_REF] Poezevara | Discovering emerging graph patterns from chemicals[END_REF], related to the data mining field, aims to discover subgraphs with a large difference of frequencies in a set of positive and negative examples. The second subfamily [START_REF] Brun | Symbolic learning vs. graph kernels: An experimental comparison in a chemical application[END_REF], more related to the machine learning field, builds a structural description of each class of molecule so that the classification is conducted by mean of a structural matching between each prototype and a graph representation of an input molecule. Both subfamilies are however mainly restricted to the classification field.

Graph kernels can be understood as symmetric graph similarity measures. Using a semi definite positive kernel, the value k(G, G ′ ) where G and G ′ encode two input graphs corresponds to a scalar product between two vectors ψ(G) and ψ(G ′ ) in an Hilbert space. Graph kernels provide thus a natural connection between structural pattern recognition and graph theory on one hand and statistical pattern recognition on the other hand. A large family of kernels is based on the definition of a bag of patterns for each graph and deduces graph similarity from the similarity between bags. Many of them are defined on linear patterns. For example, [START_REF] Kashima | Kernels for graphs[END_REF] base graph kernels definition on the comparison of sets of walks extracted from each graph. [START_REF] Ramon | Expressivity versus efficiency of graph kernels[END_REF] and [START_REF] Mahé | Graph kernels based on tree patterns for molecules[END_REF] define kernels using an infinite set of tree patterns instead of walks. These methods improve the limited expressiveness of linear features such as walks. A common drawback of kernels based on walks and tree patterns is that the global similarity between two graphs is based on an implicit enumeration of their common patterns. Such an implicit enumeration, on infinite bags, doers not allow to analyse the relevance of each pattern according to a given dataset or a given property. Using such kernels the weight of each pattern is thus set using generic functions with few parameters allowing to adapt the kernel to a given dataset. Instead of decomposing graphs into an infinite set of substructures (i.e. walks or trees), [START_REF] Shervashidze | Efficient graphlet kernels for large graph comparison[END_REF] define a kernel based on the distribution of a predefined set of unlabeled subgraphs, called graphlets. The weighting of each graphlet according to a specific dataset is also not proposed by the authors.

Another approach to the definition of graph kernels is proposed by [START_REF] Neuhaus | Bridging the gap between graph edit distance and kernel machines[END_REF] and [START_REF] Riesen | Graph embedding in vector spaces by means of prototype selection[END_REF]. This approach aims to define semi definite positive kernels from the notion of edit distance. The main challenge of this approach is that edit distance does not always fulfill all requirements of a metric and hence does not always readily lead to a semi definite positive kernel. Using non definite positive kernel, the embedding space associated to the kernel is no more an Hilbert space but a Krein space [START_REF] Ong | Learning with nonpositive kernels[END_REF] on which many usual mathematical results should be interpreted differently. Problems defined using non definite kernels usually provide a more direct formulation of the problem to minimise but should use with caution most of usual mathematical operations which must be interpreted within a Krein space. Conversely, methods based on definite positive kernel may have to regularize the kernel in order to force its definite positiveness but may apply confidently a large family of mathematical tools defined within Hilbert spaces.

Approaches based on bag of patterns and graph edit distance have complementary advantages and drawbacks. Indeed, the first approach based on the notion of bag of substructures does not readily allows to capture global information about graphs. This last point may be an advantage when a given physical property should be explained by the presence of several substructures instead of the global shape of the molecule. Conversely, properties explained by global molecule's similarities may be better captured by an approach based on graph edit distance.

This paper presents two new kernels: A first kernel, presented in Section 2, is explicitly based on the problem stated by classification or regression kernel methods. This kernel, based on both graph edit distance and graph Laplacian kernel, is regularised in order to force the definite positiveness of the kernel. These points distinguish our kernel from the ones presented, for example, by [START_REF] Neuhaus | Bridging the gap between graph edit distance and kernel machines[END_REF]. A method to update efficiently this kernel is also proposed.

Our second kernel, presented in Section 3, is related to [START_REF] Shervashidze | Efficient graphlet kernels for large graph comparison[END_REF] and [START_REF] Mahé | Graph kernels based on tree patterns for molecules[END_REF] as it is based on an enumeration of subtrees within an acyclic labeled graph. The first contribution of this kernel is an extension of [START_REF] Shervashidze | Efficient graphlet kernels for large graph comparison[END_REF] to larger and labeled subgraphs. The second contribution is, unlike [START_REF] Mahé | Graph kernels based on tree patterns for molecules[END_REF] method, an explicit enumeration of subtrees within a graph. This last point has important consequences. Indeed, while the set of potential acyclic subgraphs is very large and even infinite for infinite set of labels, our explicit enumeration of subtrees provides a finite set for each data set. This last point allows, using an appropriate variable selection method, to select the substructures which are pertinent for a given property. Such a selection step both provides useful information to interpret chemical mechanisms connected with a given property and allows to improve prediction results.

Using our approach, each molecule is described by a vector encoding the frequencies of its subtrees. From this point of view, our method may be connected with graph embedding methods. However, our vectors are embedded in a space of very large dimension. Indeed, considering a set of N differents labels, the number of different labeled trees up to size 6 is equal to Γ

1 N + Γ 2 N + Γ 3 N + 2Γ 4 N + 3Γ 5 N + 6Γ 6 N
, where Γ k N denotes the number of possible subsets with repetitions of size k among N . Since the atoms can be labeled by each of the 118 chemical elements, our feature space has a dimension higher than Γ 6 118 ≃ 4.25 × 10 9 . Nevertheless, for each dataset, we enumerate a finite set of different treelets which defines the subspace on which we encode our set of graphs. Therefore, unlike most of graph embedding methods, the final dimension of the space describing our set of molecules is determined from the data set and not a priori. The efficiency of these two approaches is finally demonstrated in Section 4 through experiments.

2 Kernel from Edit Distance

Edit Distance

An edit path between two graphs G and G ′ is defined as a sequence of operations transforming G into G ′ . Such a sequence may include vertex or edge addition, removal and relabeling. Given a cost function c(.), associated to each operation, the cost of an edit path is defined as the sum of its elementary operation's costs. The minimal cost among all edit paths transforming G into G ′ is defined as the edit distance between both graphs. A high edit distance indicates a low similarity between two graphs while a small one indicates a strong similarity.

According to [START_REF] Neuhaus | Bridging the gap between graph edit distance and kernel machines[END_REF], the computational cost of the exact edit distance grows exponentially with the size of graphs. Such a property limits the computation of exact edit distance to small graphs. To overcome this problem, [START_REF] Fankhauser | Speeding up graph edit distance computation through fast bipartite matching[END_REF] defined a method to compute an approximate edit distance in O(nd 2 ) where n and d are respectively equal to the number of nodes and to the maximal degree of both graphs.

Graph Laplacian Kernel

Unfortunately, edit distance doesn't always define a metric and trivial kernels based on edit distance may not be semi definite positive (Section 1). [START_REF] Neuhaus | Bridging the gap between graph edit distance and kernel machines[END_REF] proposed several methods to overcome this important drawback. However, the proposed kernels are not explicitly based on the minimization problem addressed by kernel methods. This minimization problem may be stated as follows: Given a kernel k and a dataset of graphs D = {G 1 , . . . , G n }, the Gram matrix K associated to D is an n × n matrix defined by K i,j = k(G i , G j ). Within the kernel framework, a classification or regression problem based on K may be stated as the minimization of the following formula on the set of real vectors of dimension n:

f * = arg min f ∈R n CLoss(f, y, K) + f t K -1 f (1)
where CLoss(., ., .) denotes a given loss function encoding the distance between vector f and the vector of known values y.

Each coordinate f i of such a vector corresponds to a value attached to graph G i . A vector f may thus be considered as a function mapping each graph of the database {G 1 , . . . , G n } to a real value. As denoted by [START_REF] Steinke | Kernels, regularization and differential equations[END_REF], the term f t K -1 f in Equation 1 may be considered as a regularization term which counter balances the fit to data term encoded by function CLoss (., ., .). Therefore, the inverse of K (or its pseudo inverse if K is not invertible) may be considered as a regularization operator on the set of vectors of dimension n and hence on the set of real functions defined on {G 1 , . . . , G n }. Conversely, the inverse (or pseudo inverse) of any semi definite positive regularization operator may be considered as a kernel. We thus follow a kernel construction scheme recently introduced [START_REF] Brun | Symbolic learning vs. graph kernels: An experimental comparison in a chemical application[END_REF] which first builds a semi definite positive regularization operator on the set of functions mapping each graph {G 1 , . . . , G n } to a real value. The inverse, or pseudo inverse of this operator defines a kernel on the set {G 1 , . . . , G n }.

In order to construct this regularization operator, let us define a n × n adjacency matrix W by:

W ij = e -d(G i ,G j ) σ (2)
where d(., .) denotes the edit distance and σ is a tuning variable. The Laplacian of W is defined as l = ∆ -W where ∆ is a diagonal matrix defined by:

∆ i,i = n j=1 W i,j (3)
Classical results from spectral graph theory [START_REF] Chung | Spectral graph theory[END_REF] establish that l is a symmetric semi definite positive matrix whose minimal eigenvalue is equal to 0. Such a matrix is thus not invertible. To overcome this problem, [START_REF] Smola | Kernels and regularization on graphs[END_REF] define the regularized Laplacian l = I + λl of W where λ is a regularization coefficient. The minimal eigenvalue of l is equal to 1 and this matrix is thus definite positive. Moreover, given any vector f , we have:

f t lf = f 2 + λ n i,j=1 W ij (f i -f j ) 2 (4)
Intuitively, minimising Equation 4 leads to build a vector f with a small norm which maps graphs with a small edit distance (and thus a strong weight) to close values. Such a constraint corresponds to the regularization term required by Equation 1 in order to smoothly interpolate the test values y over the set of graphs {G 1 , . . . , G n }. Our unnormalized kernel, is thus defined as K un = l-1 .

Note that a regularized normalized Laplacian kernel may alternatively be considered by introducing the matrix:

L = ∆ -1 2 l∆ -1 2 (5)
We have in this case, for any vector f :

f t Lf = n i=1 f 2 i ∆ ii + λ n j=1 W ij ∆ ii ∆ jj (f i -f j ) 2
The matrix L is definite positive and its associated kernel is defined as K norm = L-1 . Note that, our regularized normalized Laplacian kernel is not defined as the inverse of the regularized normalized Laplacian:

I + λ∆ -1 2 l∆ -1 2 (6)
This new formulation is however consistent with the regularization constraint which should be added to Equation 1 and provides significant advantages in the context of incoming data (Section 2.3).

Alternative regularization schemes [START_REF] Smola | Kernels and regularization on graphs[END_REF][START_REF] Fouss | An experimental investigation of graph kernels on a collaborative recommendation task[END_REF] may be applied to the Laplacian matrix. One of this regularization scheme of particular interest is the Laplacian exponential diffusion matrix d = exp(λl) associated to the kernel K = exp(-λl). This last kernel may be considered as a generalisation of the regularized Laplacian with further constraints on the derivatives of function f . However, this kernel does not allow to update efficiently kernel values in the context of incoming data (Section 2.3).

Incoming Data

Let us first consider a kernel defined from the unnormalized Laplacian. Given our learning set D = {G 1 , . . . , G n }, the test of a new graph G within a regression or classification scheme requires to update the unnormalized Laplacian l with this new graph and to compute the updated kernel defined as the inverse of the regularized and unnormalized Laplacian K = (I + λl) -1 . This direct method has a complexity equal to O((n + 1) 3 ), where n is the size of our data set. Such a method is thus computationally costly, especially for large datasets. In this section, we propose a method to reduce the complexity of this operation.

Given the regularized and unnormalized Laplacian ln = (I n +λ(∆ n -W n )) defined on the dataset D, its updated version ln+1 defined on D ∪ {G} may be expressed as follows:

ln+1 = ln -δ n B B t 1 -i B i where B = (-λexp( -d(G,G i ) σ
)) i={1,...,n} is deduced from the weights between the new input graph G and each graph (G i ) i={1,...,n} of our dataset and δ n is a diagonal matrix with (δ n ) i,i = B i .

The minimal eigenvalue of ln+1 is equal to 1 (Section 2). This matrix is thus invertible, and its inverse may be expressed using a block inversion scheme:

K un = ( ln+1 ) -1 = Γ Θ Λ Φ with          Γ = E -1 + ΦE -1 BB t E -1 Θ = -E -1 BΦ Λ = -ΦB t E -1 Φ = (1 -i B i -B t E -1 B) -1 (7
) where E = lnδ n . Note that Φ corresponds to a scalar. The computation of our new kernel, using equation 7, relies on the computation of the inverse of the matrix E = ln + δ n which may be efficiently approximated using a development to order N of (I -l-1

n δ n ) -1 : ( ln -δ n ) -1 = l-1 n (I -l-1 n δ n ) -1 ≈ N k=0 l -k-1 n δ k n (8) Such a sum converges since l-1 n δ n 2 < 1, for λ < 1. Indeed: l-1 n δ n 2 ≤ l-1 n 2 δ n 2 ≤ δ n 2 ≤ λ max i=1,n exp( -d(G, G i ) σ )
The last term of this equation is strictly lower than one for any λ lower than one. Moreover, basic matrix calculus show that the approximation error is lower than ǫ for any N greater than:

log(2ǫ) log(max i=1,n exp( -d(G,G i ) σ )) . (9)
Equation 8 allows to approximate the inverse of ( lnδ n ) by a sum of pre computed matrices l -k-1 n multiplied by diagonal matrices. Using such pre calculus, the inverse of ( lnδ n ) and hence the computation of our new kernel is achieved in N n 2 .

If we now consider the regularized normalized Laplacian (Section 2.2)

L = ∆ -1 2 l∆ -1 2 , its inverse is defined as: L-1 = ∆ 1 2 l-1 ∆
1 2 and we have:

K norm = ∆ 1 2 K un ∆ 1 2 (10)
The update of the regularized and normalized Laplacian kernel may thus be deduced from the one of the regularized unnormalized Laplacian kernel.

Treelet Kernel

Kernels based on edit distance rely on a direct comparison of each pair of graph. An alternative strategy consists to represent each graph by a bag of patterns and to deduce the similarity between two graphs from the similarity of their bags. This strategy may provide semi definite kernels hereby avoiding the necessity to regularize the whole gram matrix for each incoming data (Section 2.3). As mentioned in Section 1, most of kernels of this family are based on linear patterns (bags of paths, trails or walks). [START_REF] Shervashidze | Efficient graphlet kernels for large graph comparison[END_REF] describe a method to enumerate, for any unlabeled graph, all its subgraphs (called graphlets) composed of up to 5 nodes. We propose here to adapt this method to the enumeration of subtrees of labeled and unlabeled acyclic graphs up to size 6. The resulting patterns are called treelets (Figure 1). This bound on the treelet's size corresponds to a compromise between the expressiveness of our kernel and the time required to enumerate all types of treelets. Indeed, as shown by [START_REF] Otter | The number of trees[END_REF] the number of different unlabeled treelets grows exponentially with the number of atoms. Moreover, according to [START_REF] Isaacs | Physical Organic Chemistry[END_REF], the influence of one atom on a chemical property does not usually exceed 3 bonds. Hence, considering larger set of treelets would forbid the design of specific procedures to retrieve each type of treelet without a clear gain on the prediction of chemical properties.

G 0 G 1 G 3 G 12 G 2 G 5 G 4 G 7 G 8 G 6 G 9 G 10 G 11 G 13
Figure 1: Acyclic and unlabeled graphlets whose maximal size is equal to 6. Centers of 3-star and 4-star are surrounded

Computing Embedded Distribution

Structural Analysis

Each treelet corresponds to a subtree of some graph of our database. Such a treelet may be denoted as t = (V t , E t , µ t , ν t ) where µ t and ν t denote vertex and edge labeling functions. When dealing with unlabeled graphs, treelet enumeration consists to classify treelets according to their graph structures (V t , E t ). The cardinal of V t being bounded, only a finite set of tree structures called tree patterns may correspond to a treelet. Following [START_REF] Shervashidze | Efficient graphlet kernels for large graph comparison[END_REF], the identification of tree patterns is initiated by an enumeration of all paths whose length is lower than or equal to 6. A recursive depth first search with a maximum depth equals to 6 is thus performed from each node of the graph. Note that, using such an enumeration, each path is retrieved from its two extremities and is thus counted twice. In order to prevent this problem, each path composed of at least two nodes is counted 1 2 times. This first step provides the distribution of patterns G 0 , G 1 , G 2 , G 3 , G 4 and G 5 (Figure 1).

Our method to compute the distribution of remaining patterns is based on the detection of nodes of degree 3, 4 and 5. These nodes, respectively denoted R 3-star , R 4-star and R 5-star , are the center of 3-star, 4-star and 5star patterns. Note that a 4-star pattern (G 8 ) contains four 3-star patterns

(a) G 7 (b) Four decompositions of G 7 into G 4 Figure 2: Tree pattern G 7 contains 4 G 4 .
Pattern Source pattern Condition

G 7 3-star |{v; v ∈ N (R 3-star ); d(v) ≥ 2}| ≥ 1 G 9 3-star |{v; v ∈ N (R 3-star ); d(v) ≥ 2}| ≥ 2 G 10 3-star ∃v 0 ∈ N (R 3-star ); d(v 0 ) ≥ 2 and |{v; v ∈ N (v 0 ) -{R 3-star }; d(v) ≥ 2}| ≥ 1 G 11 4-star |{v; v ∈ N (R 4-star ); d(v) ≥ 2}| ≥ 1 G 12 3-star |{v; v ∈ N (R 3-star ); d(v) ≥ 3}| ≥ 1
Table 1: Conditions characterizing patterns derived from 3-star and 4-star. N (v) and d(v) denote respectively the set of neighbours and the degree of vertex v.

(Figure 2). This first enumeration of nodes of degrees 3, 4 and 5 provides the distribution of patterns G 6 , G 8 and G 13 . Patterns G 7 , G 9 , G 10 and G 12 are enumerated from the neighbourhood of 3-star patterns. For example, the tree pattern G 7 matches the neighborhood of a 3-star node if at least one of the neighbours of this 3-star has a degree greater than or equal to 2. Pattern G 11 is the only subtree derived from a 4-star. Properties characterizing patterns with a 3 or 4-star are summarized in Table 1. Note that pattern G 12 is symmetric since it contains two centers of 3-star. Such a treelet will thus be counted twice (once from each of its 3-star) and must be counted for 1 2 . One may easily check that no isomorphism exists between patterns depicted in Figure 1. Moreover, as shown by [START_REF] Cayley | On the analytic forms called trees, with applications to the theory of chemical combinations[END_REF], the number of different acyclic and unlabeled graphs with up to 6 vertices and a maximum degree of 4 is equal to 13. Therefore, tree patterns G 0 to G 12 (Figure 1) represent, up to isomorphisms, all unlabeled acyclic graphs whose size is lower than or equal to 6 and whose vertex degree is bounded by 4. Adding G 13 to this set provides all unlabeled acyclic graphs with a size lower than or equal to 6.

Labeling Information

Method described in Section 3.1.1 allows us to associate each treelet to a tree pattern. Using databases of unlabeled graphs, this method can be used to enumerate the number of treelets corresponding to each tree pattern. However, using databases of labeled graphs, treelets with different vertices and edge labels may be associated to a same tree pattern. We propose in this section to identify each treelet by a code composed of two parts: A structural part corresponding to the index of its tree pattern (Section 3.1.1) and a canonical key defined as a sequence of vertices and edge labels. Such a sequence is specific to each tree pattern and designed so as two treelets with a same code must be isomorphic.

Such a key is trivial for linear patterns, i.e. paths. Each path may indeed be associated to two sequences composed of alternated vertices and edge's labels encoding the two possible traversals of this path. By convention, the canonical key of such linear tree patterns is defined as the sequence with the lowest lexicographic order.

Morgan numbering

Let us consider a non linear treelet t = (V t , E t , µ t , ν t ), where µ t and ν t denote respectively vertex and edge labeling functions. Our canonical key for non linear patterns is based on the extended connectivity concept introduced by [START_REF] Morgan | The generation of a unique machine description for chemical structures-a technique developed at chemical abstracts service[END_REF]. This concept is based on an additional vertex labeling function λ from V to N called extended labeling. This function is defined by an iterative process which initialises each extended label λ(v) to the degree of v. This initial labeling is then extended by assigning to each vertex the sum of extended labels of its neighbors. This summation process is iterated while it increases the number of different labels. The resulting set of extended labels is the same for two isomorphic graphs and is unique for each tree pattern. Figure 3 shows the extended labels computed on non linear patterns.

Since two adjacent vertices v and v ′ of a treelet may be compared according to λ(v) and λ(v ′ ), our extended label defines a partial order relationship between adjacent vertices of a treelet. Following [START_REF] Morgan | The generation of a unique machine description for chemical structures-a technique developed at chemical abstracts service[END_REF], we encode this partial order relationships by a rooted tree. Treelets G 6 to G 11 have a unique local maximum of the extended labeling function. We thus root these trees on this unique vertex. Treelet G 12 has two local maxima located on vertices v 0 and v 3 (cf. Figure 3). This treelet is thus associated to two rooted trees respectively rooted on v 0 and v 3 .

v 0 v 1 v 2 v 3 v 4 v 5 G 12 1 1 1 1 3 3 v 0 v 1 v 2 v 3 v 4 v 5 G 11 4 2 1 1 1 1 v 0 v 1 v 2 v 3 v 4 v 5 G 10 2 3 5 4 3 3 v 0 v 1 v 2 v 3 G 6 3 1 1 1 v 0 v 1 v 2 v 3 v 4 G 7 1 2 3 1 1 v 0 v 1 v 2 v 3 v 4 G 8 4 1 1 1 1 v 0 v 1 v 3 v 4 v 2 v 5 G 9 5 4 4 3 2 2 G 13 v 0 v 1 v 2 v 3 v 4 v 5 1
Key Construction Scheme Our construction scheme of the canonical key of a treelet is based on a traversal of its rooted tree. The construction of our key requires to sort the children of each internal vertex of the tree in order to define a unique traversal and hence a unique key. This sorting step is achieved through the following recursive construction of our canonical key:

The key of each leaf v, key(v), is defined as the null label. For each internal node v of the tree, let us consider its set of children {v 1 , . . . , v n }. This set is first sorted according to λ(v i ) and then according to the chain defined as the concatenation of µ t (v i ), ν t (v, v i ) and key(v i ). Assuming such an order on {v 1 , . . . , v n }, the key of vertex v is defined as:

key(v) = n i=1 µ t (v i ).ν t (v, v i ) . n i=1 key(v i ) (11)
where denotes the concatenation operator.

Using such a construction scheme, the label of each vertex is encoded in the key of its father. In order to integrate the label of the root node, we define the key of a tree rooted on vertex r as µ t (r).key(r).

Treelets G 6 to G 11 are encoded by a single rooted tree, and their canonical code is defined as the index of their tree pattern concatenated with the key of their rooted tree. On the other hand, treelet G 12 is associated to two trees rooted on its two local maxima. Its canonical code is defined as the index of its tree pattern (12) concatenated with the lowest key of its two associated rooted trees according to the lexicographic order.

Note that, unlike Morgan's representation, our key does not encode the structure of the graph which is explicitly encoded by the index of the tree pattern.

Equal keys and graph isomorphism

Our canonical key is based on extended labels deduced from the structure of the treelet and from vertex and edge labeling functions. Hence, two isomorphic treelets are associated to a same canonical key. Conversely, since one may easily build a linear pattern from its canonical key, two linear patterns (G 0 to G 5 ) with a same canonical key must be isomorphic.

Within our key construction scheme, the label of a vertex with an unique extended label is located at a fixed position within the canonical key of its treelet. The label of such a vertex may thus be retrieved without any ambiguity from a canonical key. However, a set of children {v 1 , . . . , v n } of a same parent node v with a same extended label will be sorted according to the sequence of edge and vertex labels µ t (v i )ν t (v, v i )key(v i ). This sorting step provides a unique label for two isomorphic treelets but does not allow us to distinguish between any permutations of vertices {v 1 , . . . , v n }. We have thus to check, for each treelet, that all permutations of vertices allowed by our code correspond to isomorphic treelets.

Vertices with a same extended label within tree patterns G 6 , G 7 , G 8 , G 10 , G 11 and G 12 are represented by black squares ( ) in Figure 3. For each tree pattern, using our key construction scheme, these vertices of degree one are the child of the unique vertex they are connected to. Therefore, our key does not distinguish any permutation among these vertices. Note that, since these vertices have a degree one and are connected to a same vertex, any permutation swapping two of these vertices leads to an isomorphic treelet.

The rooted tree associated to treelet G 9 does not allow us to distinguish between branches v 2 v 3 and v 4 v 5 (Figure 5(b)) since both v 2 and v 4 have an extended label equal to 4. However, the simultaneous permutation of v 2 with v 4 and v 3 with v 5 provides an isomorphic graph. In the same way, the canonical key of G 12 (Figure 5(c)) does not allow us to distinguish permutations between v 0 and v 3 nor permutations between v 1 , v 2 on one hand and v 4 , v 5 on the other hand. However, as illustrated in Figure 5(c), all these permutations provide isomorphic treelets. Treelets G 0 to G 12 are thus uniquely identified by their codes up to isomorphisms.

v 0 v 1 v 3 v 4 v 2 v 5 5 4 4 2 2 3 (a) G 9 v 0 v 1 v 4 v 2 v 3 v 5 (b) Rooted tree of G 9 1 1 3 3 1 1 v 0 v 1 v 2 v 3 v 4 v 5 (c) G 12 v 2 v 5 v 4 v 0 v 1 v 3 (d) First rooted tree of G 12 v 2 v 5 v 4 v 0 v 1 v 3 (e) Second rooted tree of G 12

Complexity

Our treelet enumeration being restricted to tree patterns, its complexity for the enumeration of treelets up to size 5 is bounded by the one of [START_REF] Shervashidze | Efficient graphlet kernels for large graph comparison[END_REF] : O(nd 4 ), where n is the number of nodes and d is the maximum degree of the graph. Using a depth first search, unlabeled linear patterns can be enumerated in O(nd k-1 ), the enumeration of paths of length 6 may thus be achieved in O(nd 5 ). The detection of p-stars requires to enumerate all subsets of p neighbours of each node of the graph. This enumeration is achieved in O(nd p ) operations. The enumeration of G 6 , G 8 and G 13 is thus respectively performed in O(nd 3 ), O(nd 4 ) and O(nd 5 ) operations. The enumeration of and the enumeration of all 4-stars is performed in O(nd 4 ), the enumeration of all treelets corresponding to G 11 is performed in O(nd 5 ). Therefore, the enumeration of all unlabeled treelets is performed in O(nd 5 ). Since each unlabeled treelet is finite, the complexity required to compute its key is performed in O(1). Therefore, the overall complexity required to enumerate labeled treelets from a graph remains equals to O(nd 5 ). Note that considering bounded degree graphs, the complexity required to enumerate treelets is linear with the number of graph's nodes.

v 4 v 2 v 0 v 1 v 5 v 3 (a) G 9 v 0 v 1 v 2 v 3 v 4 v 5 (b) G 10 v 3 v 0 v 1 v 2 v 4 v 5 (c) G 12 v 0 v 1 v 2 v 3 v 4 v 5 (d) G 11

Treelet Kernel Definition

Based on our treelet enumeration, we define a function f which associates to each graph G a vector f (G) encoding the distribution of its treelets. Each component of this vector, denoted the spectrum of G, is equal to the frequency of a given treelet t in G:

f (G) = (f t (G)) t∈T (G) with f t (G) = |(t ⊆ G)| (12)
where T (G) denotes the set of treelets of G.

Given the function f , the kernel between two graphs is defined as a sum of kernels between the different number of treelets common to both graphs:

k T reelet (G, G ′ )q = t∈T (G)∩T (G ′ ) K (f t (G), f t (G ′ )) ( 13
)
where T (G) ∩ T (G ′ ) denotes the set of common treelets between G and G ′ . The kernel K between f t (G) and f t (G ′ ) corresponds to a definite positive kernel between real numbers. In our experiments (Section 4), the RBF, scalar and polynomial kernels have been tested in order to select the best one on each dataset.

Treelet Weighting

Among the set of treelets found within a dataset, some of them may be irrelevant to explain a given property. Such irrelevant treelets may induce useless calculus and decrease performances of regression or classification algorithms.

A first trivial approach to this problem considers all possible permutations over the set of treelets. This approach implies to test 2 p different sets of treelets, where p is the number of different treelets of the dataset. Such a brute force approach is intractable even for small values of p.

In order to define a set of relevant treelet in a less greedy way, we propose two iterative approaches defined for regressions problems. The first one, called forward selection, is an iterative method which starts with an empty set and adds one treelet at each step. The chosen treelet is the one which produces the best regression result (Alg. 1) on the training set using a leave one out procedure. To evaluate the quality of the result, we use the Residual Sum of Squares (RSS) defined as the sum of squared prediction error made for each molecule. A similar approach consists to start from the whole set of treelets and to remove one treelet at each iteration. This second approach is called backward elimination (Alg. 2). Both methods imply to test p(p+1) 2 set of treelets. The main difference between these two approaches concerns groups of treelets which are strongly related to a given property when considered simultaneously. For such groups, the removal of any treelet by backward elimination induces a strong increase of the prediction error. These groups are thus preserved by backward elimination algorithm. On the other hand, such groups may not be created by forward selection since the addition of a single treelet may not induce a significant decrease of the prediction error.

Experiments

Classification Problem

Our first experiment evaluates graph kernels defined in Section 2 and 3 on a classification problem. This problem is defined on the monoamine oxidase (MAO) dataset1 which is composed of 68 molecules divided into two classes: 38 molecules inhibit the monoamine oxidase (antidepressant drugs) and 30 do not. These molecules are composed of different chemical elements and are thus encoded as labeled graphs. Classification accuracy is measured for each method using a leave one out procedure with a two-class SVM. This classification scheme is made for each of the 68 molecules of the dataset.

Table 2 shows results obtained by different methods. The first two lines correspond to methods based on linear patterns. [START_REF] Suard | Kernel on bag of paths for measuring similarity of shapes[END_REF] (Line 1) computes the similarity between two graphs from the average similarity be-

Method Classification Accuracy

(1) [START_REF] Suard | Kernel on bag of paths for measuring similarity of shapes[END_REF] 80% (55/68) (2) [START_REF] Vishwanathan | Graph Kernels[END_REF] 82% (56/68) (3) [START_REF] Neuhaus | Bridging the gap between graph edit distance and kernel machines[END_REF] 90% (61/68) (4) [START_REF] Riesen | Graph embedding in vector spaces by means of prototype selection[END_REF] 91% (62/68)

( tween each pair of paths extracted from both graphs. Vishwanathan et al.

(2010) (Line 2) counts the number of identical random walks of two graphs. Line 3 [START_REF] Neuhaus | Bridging the gap between graph edit distance and kernel machines[END_REF] corresponds to a Gaussian kernel on the suboptimal edit distance and Line 4 [START_REF] Riesen | Graph embedding in vector spaces by means of prototype selection[END_REF] corresponds to an embedding of the graphs into a vector encoding the edit distance between the graphs and a set of selected prototypes. Lines 5 and 6 correspond to graph Laplacian kernel (Eq. 10) using a suboptimal graph edit distance [START_REF] Fankhauser | Speeding up graph edit distance computation through fast bipartite matching[END_REF] with node substitution and edge deletion costs set to 1 and edge substitution cost set to the sum of incident node substitution costs. Then, the last two lines correspond to methods based on non linear patterns. [START_REF] Mahé | Graph kernels based on tree patterns for molecules[END_REF] (Line 7) deduce the similarity between two graphs from the number of common tree-patterns and treelet kernel (Line 8) is defined in Section 3. In this experiment, this last kernel is defined using a RBF kernel between treelet's frequencies (Equation 13).

Graph Laplacian kernel methods obtain a classification accuracy of 90%. We can note that others prediction models based on edit distance obtains similar results: (Riesen et al., 2007) obtains 91% and(Neuhaus andBunke, 2007)'s obtains 90% of classification accuracy. This last kernel may however be non semi definite positive. We may additionally notice that the use of our fast inversion method (Table 2, Line 6) defined in Section 2.3 does not

Method Classification Accuracy

(1) [START_REF] Riesen | Iam graph database repository for graph based pattern recognition and machine learning[END_REF] 97.3% (2) [START_REF] Suard | Kernel on bag of paths for measuring similarity of shapes[END_REF] 98.5% (3) [START_REF] Vishwanathan | Graph Kernels[END_REF] 98.5% (4) [START_REF] Neuhaus | Bridging the gap between graph edit distance and kernel machines[END_REF] 99.7% (5) [START_REF] Riesen | Graph embedding in vector spaces by means of prototype selection[END_REF] 98.2% (6) Graph Laplacian Kernel 99.3% (7) Treelet Kernel 99.1% 2, Line 5). The number of iterations required by this fast inversion method is determined by Equation 9. Our experiments performed on the MAO database show that a value of ǫ equals to 10 -4 induces a maximum of 9 iterations hence allowing to update the gram matrix in O(9n 2 ) instead of O(n 3 ) using a standard matrix inversion method. The low value of n on this dataset (n = 68) does not induce an important gain on execution time since the average time to update a Gram matrix using method described in Section 2.3 is 0.273ms on the MAO database while this time is equal to 0.498ms using a direct matrix inversion. The ratio between both execution times is nevertheless about 1.8 hence showing a significant gain. Thanks to the fact that they take into consideration more structural information than methods based on linear patterns, Mahé and Vert's kernel (Line 7) and our treelet kernel (Line 8) obtain the best results. Note that our method enumerates 153 differents treelets. Mahé and Vert kernel's best results are obtained when considering large and complex trees, i.e. when setting maximum depth and penalization parameters to high values.

Our second classification experiment has been performed on a graph database provided by [START_REF] Riesen | Iam graph database repository for graph based pattern recognition and machine learning[END_REF]. This database defined from the AIDS Antiviral Screen Database of Active Compounds is composed of 2000 chemical compounds some of them being disconnected. These chemical compounds have been screened as active or inactive against HIV and they are split into three different sets:

• A train set composed of 250 compounds used to train SVM.

• A validation set composed of 250 compounds used to find parameters giving the best accuracy result.

• A test set composed of remaining 1500 compounds used to test the classification model.

Table 3 shows results obtained by different methods on this dataset. Note that results obtained by [START_REF] Mahé | Graph kernels based on tree patterns for molecules[END_REF] are not displayed for this dataset since the source code provided by the authors is restricted to molecules with a degree bounded by 4. The method presented by Riesen and Bunke (Table 3, Line 1) uses the graph edit distance and k-Nearest Neighbour classifier. Others methods based on edit distance, Neuhaus and Bunke (2007) (Table 3, Line 4), [START_REF] Riesen | Graph embedding in vector spaces by means of prototype selection[END_REF] (Table 3, Line 5) and graph Laplacian kernel (Table 3, Line 6) improve accuracy of the prediction model and obtain better results than simple graph edit distance. Gaussian kernel on edit distance (Table 3, Line 4) combined with SVM obtains the best results on this dataset. Note that the regularization added to Gaussian edit distance kernel affects prediction accuracy while keeping good results. However, fast Gram matrix update hasn't be tested in this dataset since best results are obtained using a parameter λ > 1 (Section 2.3). Then, Treelet Kernel (Table 3, Line 7), based on non linear patterns and a scalar product kernel between frequencies, outperforms methods based on linear patterns (Table 3, Lines 2 and 3) but without reaching the accuracy of methods based on edit distance. Note however, that 4875 different treelets are enumerated from the train set.

A treelet selection procedure for classification applied on this dataset should improve prediction accuracy.

Regression Problem

Our first regression experiment is based on a database of alkanes (Cherqaoui and Villemin, 1994). An alkane is an acyclic molecule solely composed of carbons and hydrogens. A common encoding consists to implicitly encode hydrogen atoms using the valency of carbon atoms. Such an encoding scheme allows to represent alkanes as acyclic unlabeled graphs. The alkane dataset is composed of 150 molecules associated to their respective boiling points.

Using the same protocol than Cherqaoui and Villemin (1994), we evaluate the boiling point of each alkane using several test sets composed of 10% of the database, the remaining 90% being used as training set.

Method

Average Error RMSE ( • C)

( • C)
(1) Cherqaoui and Villemin (1994) 3.11 3.70 (2) [START_REF] Neuhaus | Bridging the gap between graph edit distance and kernel machines[END_REF] 5.42 10.01 (3) [START_REF] Riesen | Graph embedding in vector spaces by means of prototype selection[END_REF] 5.27 7.10 (4) [START_REF] Suard | Kernel on bag of paths for measuring similarity of shapes[END_REF] 4.66 6.21 (5) [START_REF] Vishwanathan | Graph Kernels[END_REF] 10 The first line of Table 4 shows results obtained by a method based on neural networks. The other lines correspond to methods based on graph kernels and boiling point prediction was performed using a kernel ridge regression. Poor results obtained by methods based on edit distance, including graph Laplacian kernel, can be explained by the lack of local information when dealing with unlabeled graphs. Indeed, using such graphs, the heuristic used to approximate graph edit distance [START_REF] Fankhauser | Speeding up graph edit distance computation through fast bipartite matching[END_REF] maps the set of vertices of both graphs using uniquely the degree of vertices. Such a method thus considers several mappings as equivalent if several vertices with a same degree exist in both graphs. In this case, the suboptimal graph edit distance thus induces a poor graph discrimination. This lack of local information within unlabeled graphs also explains the poor results obtained by [START_REF] Suard | Kernel on bag of paths for measuring similarity of shapes[END_REF] and [START_REF] Vishwanathan | Graph Kernels[END_REF]. Indeed, using unlabeled graphs, these kernels are based on linear structures which are only discriminated by their lengths. On the other hand, kernels based on non linear patterns, treelet kernel and [START_REF] Mahé | Graph kernels based on tree patterns for molecules[END_REF], outperform previous results of Cherqaoui and Villemin (1994) based on neural networks combined with chemical descriptors. Note that we didn't select relevant patterns in this experiment since our set of unlabeled treelet is composed of only 13 treelets. In this experiment, the kernel between treelet's frequencies, on which is based our treelet kernel, is set to the RBF kernel.

Method

Average Error ( Cherqaoui et al. (1994a) 3.01 5.102 (2) [START_REF] Neuhaus | Bridging the gap between graph edit distance and kernel machines[END_REF] 6.84 9.04 (3) [START_REF] Riesen | Graph embedding in vector spaces by means of prototype selection[END_REF] 7.39 9.94 (4) GLK (Diffusion Process) 7.27 9.84 (5) [START_REF] Suard | Kernel on bag of paths for measuring similarity of shapes[END_REF] 7.28 12.37 (6) [START_REF] Vishwanathan | Graph Kernels[END_REF] 14.64 18.95 (7) [START_REF] Mahé | Graph kernels based on tree patterns for molecules[END_REF] 5 The second regression problem is based on a dataset composed of 185 acyclic molecules (Cherqaoui et al., 1994a). This problem also consists in predicting boiling point of molecules. In contrast with the previous dataset, these molecules contain several hetero atoms and are thus represented as acyclic labeled graphs. In order to compare our results, we first use a leave one out method to predict boiling points, as described in Cherqaoui et al. (1994a). Our variable selection step requires O(p 2 ) regressions, where p denotes the number of different treelets enumerated on a dataset. This step requires about 4 hours on this dataset using a standard desktop computer. Note that this step is performed during the training set and does not influence the time required to predict the property of a molecule. However, due to this important computational time, we performed in this experiment a single variable selection step on the whole dataset of 185 molecules, instead of considering all subsets of 184 molecules as required by the leave one out protocol.

• C) RMSE ( • C) (1)
The first line of Table 5 shows results obtained using a neural network based on the embedding frequencies of 20 substructures chosen by a chemical expert (Cherqaoui et al., 1994a). Then, the next three lines show results obtained by methods based on the notion of edit distance. The Gram matrix associated to the Gaussian kernel defined on suboptimal edit distance (Line 2) obtains intermediate results. Line 3 shows results obtained by graph embed-ding method using graph edit distance and prototype selection as defined in [START_REF] Riesen | Graph embedding in vector spaces by means of prototype selection[END_REF]. The graph Laplacian kernel (Line 4) with diffusion process regularization (K = e -λl , Section 2.2) obtains similar results than alternative methods based on edit distance. We should however note that the regularization of the Laplacian induces a slight decrease of classification results obtained by [START_REF] Neuhaus | Bridging the gap between graph edit distance and kernel machines[END_REF] on this dataset. The performances of the regularized Laplacian (K = (I + λl) -1 , Section 2.2) on this dataset are below the one of the graph Laplacian kernel with diffusion process regularization and are thus not displayed in Table 5. Methods based on linear patterns [START_REF] Suard | Kernel on bag of paths for measuring similarity of shapes[END_REF][START_REF] Vishwanathan | Graph Kernels[END_REF] (Lines 5 and 6) suffer from the lack of expressiveness of linear patterns and obtain poor results when predicting the boiling point of such molecules.

Conversely, methods based on non linear patterns, Mahé and Vert (2008) (Line 7) and treelet kernel using a RBF kernel between treelet's frequencies (Line 8), outperform methods based on linear patterns but without reaching the efficiency of Cherqaoui et al. (1994a). This lack of efficiency can be explained by the number of different treelets induced by this dataset. Indeed, treelet kernel method enumerates 142 different treelets in this dataset and some of them are not related with the property to predict. Therefore, information brought by these treelets can be assimilated to noise which decreases prediction results. These prediction results can be improved by applying one of the two treelet selection methods described in Section 3.3. As shown in Table 5, prediction errors made using these selection methods (Lines 9 and 10) outperform results obtained by alternative regression methods on this dataset. Note that forward selection method reduces the set of features to 26 treelets and backward elimination reduces it to 56 treelets. Backward elimination obtains the best result by taking into account information brought by combinations of treelets (Section 3.3). An important difference between treelet kernel approach and [START_REF] Mahé | Graph kernels based on tree patterns for molecules[END_REF] is that treelet kernel method computes explicitly the distribution of each treelet within a molecule. This explicit enumeration allows us to weight each treelet separately while [START_REF] Mahé | Graph kernels based on tree patterns for molecules[END_REF] only weight features according to their depth or branching cardinality. This difference explains the better results obtained by our treelet kernel combined with a backward selection.

Table 6 shows results obtained on the same dataset using 90% of the dataset as train set and remaining 10% as test set. The method of Cherqaoui et al. (1994a) is not displayed in this table since only results based on a leave one out procedure are available. On this experiment, our treelet kernel is based on a RBF kernel between treelet's frequencies and our variable selection step is performed only using the train set. With this restricted training set, kernels based on non linear patterns still outperform kernel based on linear ones and our treelet kernel combined with a backward selection step obtain the lowest prediction error.

The first column of Table 7 shows computation times required to compute the Gram matrices of the different kernels studied in this section on our dataset of labeled acyclic molecules (see Tables 5 and6 to evaluate the performances of these kernels). The second column of Table 7 describes the mean time required to predict the boiling point of each molecule on the same data set using a leave one out procedure. Thanks to the low bounded degree of molecular graphs, treelet enumeration can be performed efficiently (Table 7, Line 7). In addition, treelet enumeration is performed once for each graph and only a sum of RBF kernels between treelet's frequencies is performed for each pair of graphs. Conversely, using an implicit enumeration, we have to perform two enumerations for each pair of graph. Note that the mean prediction time required by Graph Laplacian Kernel is similar to the one required to compute the whole Gram matrix since we have to inverse the regularized Laplacian matrix to compute the updated kernel.

Conclusion

We have proposed two new graph kernels using different approaches. The first one is based on the notion of edit distance and is adapted to cases where the similarity between molecules must consider the global shape of the molecule. Thanks to different regularization methods, we define a valid graph kernel from a suboptimal edit distance and we propose a solution to improve the computational complexity of this kernel in case of incoming data for one type of regularisation. The results obtained by this kernel are equal or slightly below the one obtained by alternative methods based on the edit distance which do not perform a regularization step. This results confirms the fact that a regularization step slightly transforms the problem to be minimized by kernel methods. However, such a regularization step is the price to pay in order to work in an Hilbert space and have a clear interpretation of all the mathematical operations performed by kernel methods. Note, nevertheless that some kernel methods have been specifically designed to work on Krein space with a clear interpretation on their limited application domain.

Our second kernel is based on the decomposition of a graph into a set of distinct substructures called treelets. Such a kernel is more specifically designed to capture the similarity between molecules sharing many sub structures. In contrast with many existing graph kernel methods only based on linear patterns, our method is based on linear and non linear patterns which provide additional local information. Moreover, our explicit enumeration of treelets allows both to define an efficient kernel since we avoid an explicit comparison of all pairs of treelets of two graphs and it additionally allows us to integrate within our kernel's definition a treelet selection step which improves the accuracy of our regression results. Our treelet kernel indeed outperforms all other kernels on all our regression experiments. This treelet selection step is however, not yet defined for classification methods and our kernel obtains on these kind of experiments good results which are however always outperformed by an alternative kernel. Designing a specific treelet selection step for classification tasks is one of the priority of our future work. Other perspectives include the extension of our treelet kernel in order to encode ring information while keeping an efficient enumeration.

Figure 3 :

 3 Figure 3: Patterns with Morgan's extended labeling. Possible permutations are represented by square nodes. Different permutations within a same treelet are encoded by different gray levels.

Figure 4 :

 4 Figure 4: Rooted trees associated to treelets based on the extended labeling function λ. Values of λ is indicated next to each vertex in (a) and (c).

Figure 5 :

 5 Figure 5: Neighbourhood analysis required to enumerate treelets of size 6 based on 3 and 4-stars.

Table 2 :

 2 Classification accuracy on the monoamine oxidase (MAO) dataset.

	5)	Normalized Standard Graph Laplacian Kernel	90% (61/68)
	(6)	Normalized Fast Graph Laplacian Kernel	90% (61/68)
	(7) Mahé and Vert (2008)	96% (65/68)
	(8) Treelet Kernel	94% (64/68)

Table 3 :

 3 Results on AIDS dataset.

modify graph Laplacian kernel's classification accuracy (Table

Table 4 :

 4 Boiling point prediction on alkane dataset. RMSE stands for Root Mean Squared Error.

	.61	16.28

Table 5 :

 5 Boiling point prediction on acyclic molecule dataset. TK stands for treelet kernel; GLK for Graph Laplacian Kernel.

	.53	9.50

Table 7 :

 7 Empirical execution times measured on the acyclic molecule regression problem.

Except AIDS dataset, all databases in this section are available on the TC15 Web page: http://www.greyc.ensicaen.fr/iapr-tc15/links.html#chemistry


B Farbos

J.-M Leyssale

J.-P Da Costa

P Weisbecker

G L Vignoles

Atomistic models of pyrolytic carbons obtained with the IGAR method

Introduction

Pyrolytic carbons (pyCs) matrices are, with carbon fibers, one of the two components in high performances carbon/carbon thermostructural composites. PyCs are turbostratic carbon materials made of a stacking of faulted graphene layers. In a recent publication [START_REF] Leyssale | Atomistic models of pyrocarbons: reconstruction procedure and structural features[END_REF], Leyssale et al. have studied the atomic structure of a Rough Laminar (RL) pyCs using the Image Guided Atomistic Reconstruction (IGAR) method and showed that in addition to a high percentage of nanosized disoriented graphene domain, defined by threefold atoms belonging only to hexagonal rings (called pure C 6 atoms), defects are mainly made of intrasheet grain boundaries made of lines or networks of pentagon/heptagon pairs and of screw dislocations connecting different sheets of a stack to each other through in a kind of "car park access ramp" arrangement. A comparison of AP and HT (inducing a partial graphitization) models shows that HT model have smoother fringes (related to a higher sp 2 fraction), larger monolayer aromatic domains, larger C 5 /C 7 networks and fewer dislocations. These results are sum up in the table 1 Very recently, Weisbecker et al. [START_REF] Proceeding | [END_REF] have characterized two kinds of pyCs. A rough Laminar (RL) pyC, that are pyCs with a high anisotropy (R A = 8.2±0.5) and with few defects (FWHM D = 83±2 cm -1 ), and a regenerated Laminar (ReL) pyC, that are pyCs with similar anisotropy (R A = 7.0±1.0) as RL but with more defects (FWHM D = 204±20 cm - 1 ). Obviously, the differences in texture and structure influence their properties [START_REF] Vallerot | M atrice de pyrocabone : propriétés, structure et anisotropie optique[END_REF]. This study aims at obtaining atomistic representations of these RL and ReL pyCs and at describing their nano-texture/structure. We will first describe some improvements in the IGAR method. Then, the nature of defects in these materials, in terms of non-hexagonal rings and screw dislocations, as well as their hydrogen atoms contents and locations, will be thoroughly compared and discussed. These models will also be validated against experimental HRTEM images and pair distribution functions (PDF) (and structure factors) obtained from neutron diffraction.

Methods

As explained by Leyssale et al. [START_REF] Leyssale | Atomistic models of pyrocarbons: reconstruction procedure and structural features[END_REF][START_REF] Leyssale | An image guided atomistic reconstruction of pyrolitic carbons[END_REF], atomistic models of laminar pyCs are constructed using the combination of 2D HRTEM lattice fringe image analysis, 3D image synthesis and atomistic simulated annealing. Using a parametric approach [START_REF] Costa | Synthesis of solid textures based on a 2D example : application to the synthesis of 3D carbon structures observed by transmission electronic microscopy[END_REF], a synthetic HRTEM-like 3D image is obtained from an anisotropic and statistically identical HRTEM 2D experimental image. The IGAR me thod consists in a quench molecular dynamics simulation during which atoms are guided by two potentials. One is the second generation reactive empirical bond order potential U REBO based on short range interatomic interactions describing C-C bonds. The other is an external potential U HRTEM based on the synthetic HRTEMlike 3D image leading the atoms toward the fringes of the image. After this quench, the system is heated to 300K, and a simulation controlled by an U AIREBO potential, derived from U REBO potential including an adaptive treatment of van der Waals interactions between non-bonded atoms, is carried out, without the constraint of the image, to allow the relaxation of the structure between the graphene sheets.

To obtain both structurally more converged and larger models in a reasonable calculation time, our IGAR computer code was first improved using the Cell linked-lists method (consisting in dividing the simulation box in many cells, so neighbors of a given atom are only calculated in its cell and its neighboring cell rather than in the entire box, during the update of the Verlet lists) [START_REF] Allen | Computer Simulation of Liquids[END_REF] and the parallel programming on shared memory systems (via openMP). The atomistic construction in this method is costly probably because of the very long liquid quench simulations permitting the reorganization of the atoms. So, the ramp of the quench rates was also improved by minimizing the U REBO potential during the quench.

Results and discussion

Thanks to our improvements, the calculation time was be divided by about 20 and it now increases almost linearly with the number of atoms. Preliminaries studies were carried out on RL pyCs as prepared (AP) and after a heat treatment (HT) previously studied in [START_REF] Leyssale | Atomistic models of pyrocarbons: reconstruction procedure and structural features[END_REF]. As seen in fig. 1, the ramp was improved in determining temperature where the potential decreases drastically in varying the quench rate. The structral analysis corresponding to these ramps is reported on table 2.

Table 2. Analyse structurale trempe AP et HT

Structural analysis of these new models shows that the structural convergence is reached for the carbons having 2, 3 and 4 neighbors (named respectively « Csp, Csp 2 and Csp 3 ») and that, for example for the pyC AP using our best ramp, the Csp 2 (97.6%) are much more than Csp 3 (2.0%) and Csp (0.4%). However, the percentages of atoms belonging to rings having 5, 6, 7 and 8 carbons among the Csp 2 (named respectively « C 5 , C 6 , C 7 , C 8 ») haven't converged except C 8 (0.1%), even if the percentage of C 6 (95.3%) is close to the theoretical maximum. As already noticed in a previous study [START_REF] Leyssale | Atomistic models of pyrocarbons: reconstruction procedure and structural features[END_REF], the percentages of C 5 (2.2%) and C 7 (2.5%) are similar for every ramp, permitting to suppose the presence of pairs of defects. The percentages of pure C 6 atoms (87.7%) and pure C 6 rings (85.2%) are also not converged but high. For the pyCs HT, the percentage of Csp 2 is higher (99.0%) than the corresponding AP (so it seems Csp 2 is related to the flatness of the sheets), and there is the same quantity of C 6 , but more pure C 6 atoms and pure C 6 rings. 3D models permit to visualize the distribution of the defects (i.e. the non pure C 6 atoms) and of the graphene domains (i.e. the pure C 6 atoms). Most C 5 and C 7 are not found as isolated defects but usually merge as C 5 /C 7 pairs, as well as lines or networks of C 5 /C 7 pairs, defining some kind of intrasheet grain boundaries between nanosized disoriented aromatic domains. Then the Csp 3 defects make links between the fringes and the Csp defects locate at the edge of the domains. We can also see that different sheets of a stack are connected to each other through many screw dislocations in a kind of "car park access ramp" arrangement as already remarked [START_REF] Leyssale | Atomistic models of pyrocarbons: reconstruction procedure and structural features[END_REF]. A comparison of the AP and HT models shows that HT models have smoother fringes, larger monolayer aromatic domains, larger C 5 /C 7 networks and fewer dislocations. A quantitative analysis of the distribution of the fragments size shows some small fragments and principally one big fragment for AP whereas HT has intermediate fragments sizes for HT, proving that the density of screw dislocation is higher for AP compare to HT. Then a quantitative analysis of the average number of atoms by ring in the fragments shows that this number is higher for AP proving that the fragments shape is more threadlike for AP compare to HT. A visual analysis of the distribution of the C-C bonds lengths shows that the faulted areas (except the Csp) are also the most curved (because they correspond to the longer C-C bonds lengths). So defects create local disturbances in the fringes. But there is also curved areas without any defect, before et after relaxation of the structure. However, after relaxation of the structure (using either an AIREBO or a REBO potential) the curvatures corresponding to the defects stay but the others seem to be better distributed i.e. there is smaller curvatures in more areas rather than few big curvatures. So, the distribution of the bond lengths permits the identification of the areas under constraints. The identification of the peaks of the calculated PDF and distribution angles permits to confirm some differences between the AP and HT pyCs. Indeed, it can be easily seen that HT has more C 6 (sp 2 ), the Csp 3 are more in AP than HT, and the Csp are few. Parameters (scale number, orientation number, neighboring pixels number, and iteration number) of the images HRTEM-3D synthesis algorithm have been optimized by comparing statistically (due to the 3D synthesis) resulting synthetic 3D images with the input experimental images , and both images show a good agreement. After these preliminaries steps, we have studied the new RL pyC, named pyC1, as prepared (at aboutT=1000°C during 300h). For this, we have synthetized a 3D image of this pyC1 of 6.2 nm sides i.e. about 1.5 times L a10 (3.9±0.3 nm) and L c (4.2±0.1 nm) (because the model must be representative). The 2D filtered experimental (fig. 2 a)) and synthetic 3D images (fig. 2 b)) show a good agreement. The density d (2.12±0.01 g/cm 3 ) measured from He pycnometry permits to specify a number of atomes in the simulation box of 25273 for a pure carbon model.

The structural analysis of the obtained model (fig. 2 c)) shows results rather similar to the previous pyC AP, in terms of sp, sp 2 and sp 3 fractions (table 3), but also in terms of density of dislocation and fragment shapes … Table 3. Structurale analyse of pyC1.

The reduced experimental PDF (obtained from neutron diffraction) and calculated (from the model) until 30Å show an excellent agreement [START_REF] Proceeding | [END_REF]. However, some differences exist. Indeed, peak locations are a little shifted (actually, scaled) to higher distances when comparing model PDF to the experimental one. This "r scaling" can be explained by the fact that the model is not at mechanical equilibrium, but contains residual stresses, concomitantly due to the imposed nanotexture (the 3D HRTEM-like image) and to the small sizes of our systems (not allowing enough freedom for relaxation). Also, the peaks intensity is higher for the model. This can be explained by the fact that the model is more organized and because of residual stresses reduce freedom of bonds vibrations. We note the absence of the peak corresponding to the first interfringe distance d 002 showing the turbostratism of the pyC. The structure factors permit in particular to have acces to these interfringe distances (d 002 , …). Experimental and calculated structure factors show a good agreement. However, we note in particular that the 10 band (whose the FWHM permit to me asure L a10 ) is more intense for the model so its fringes are probably smoother. The statistical comparison between the filtered experimental and the simulated HRTEM images shows a rather good agreement, except that fringes are smoother for the model (fig. 2 a) and d)). So, in the aim to improve the model, larger models must be used and a way to quench and/or relax the structure without stretching it must be found. That's why simulated image was also obtained after a relaxation using a REBO potential, but this image seems to give similar image that with the AIREBO potential. We note that despite of these differences, we can clearly distinguish the differences between experimental or simulated images of two differents pyCs like for example pyCR1 and pyCR2 [START_REF] Leyssale | Atomistic models of pyrocarbons: reconstruction procedure and structural features[END_REF]. According to ERDA and SIMS measurements, pyC1 contains 0.7±0.4 atomic percent of hydrogen [START_REF] Weisbecker | M icrostructure of pyrocarbons from pair distribution function analysis using neutron diffraction[END_REF]. But this is in contradiction with former results (up to 6 % at.) on similar carbons, as obtained by ERDA. Adding 5% and 0.7% atomic percent of hydrogen in the models, we notice that hydrogen in graphene position (in the plan of the layers, by opposition to the ones in graphane position) become the most numerous (that is a more realistic model) on one hand when we use an hydrogen percentage of 0.7% rather than 5% (64% vs 43% H graphene) and on the other hand, when we force hydrogen to locate in the fringes thanks to the HRTEM image potential (56% vs 43% H graphene). On the 3D models, we see that hydrogen in graphene position locate at the edges of the domains (favourable) whereas the hydrogen in graphane position locate between fringes (unfavourable).

On the same way, we have studied a ReL pyC as prepared (in the range 1000-1100°C during 10h) named pyC3 (pyC2 will be a LL pyC [START_REF] Weisbecker | M icrostructure of pyrocarbons from pair distribution function analysis using neutron diffraction[END_REF]). We have synthetized a 3D image of pyc3 of 6.2nm of sides i.e. about 2 times L a10 (2.8±0.3 nm) and L c (2.6±0.1 nm). The 2D filtered experimental and synthetic 3D images show a good agreement. The density d (2.11±0.01 g/cm 3 ) permits to specify a number of atomes in the simulation box of 25154 for a pure carbon model. HRTEM images show that like pyC1, pyC3 has large anisotropic domains and undulated (FWHM D is related to the undulations of the layers associated to the low energy defects (before 2200°C)) and highly faulted 002 fringes. The anisotropy is similar for pyC1 and pyC3 but the Fourier Transform shows that crystallinity seems to be better for pyC1.

(L a10 and L c are high energy (after 2200°C) structural parameters obtained from neutron diffraction corresponding respectively to the coherence lengths along the graphene sheets and in the stacking direction. In particular, L a10 can be easily related to the size of the graphene domains). So pyC1 has larger domains and smoother fringes. Structural analysis, PDF, S(q), simulated image. … Comparison RL and ReL pyCs … à completer Table 4. Structurale analyse of pyC3.

Conclusions

The IGAR algorithm was optimized and parallelized and the ramp of the quench rates was improved. So, larger and more structural converged atomistic models of laminar (RL and ReL) pyrolytic carbons were obtained. Texture/structure of models was described by characterizing many kinds of defects and graphene domains (i.e. carbons having 2, 3 and 4 neighbors, atoms in rings of 5, 6, 7 and 8 carbons, pure C 6 atoms and pure C 6 rings, size and shape of the fragments, screw dislocations, distribution of the bonds lengths, hydrogen). Some of these structural factors have converged and others seem to be close to. The most numerous of H locates favourably if its percentage is well chosen and if the HRTEM image potential is applied on hydrogen. Models are validated both by experimental and simulated PDF, structure factors, and HRTEM images. More, the few differences permit to find futures improvements in IGAR process. In particular, larger models will be used and a way to quench and/or relax the structure without stretching the models will be found. Large and accurate models of rough laminar and regenerated laminar pyCs were obtained, studied and compared … à compléter Next, ReL pyC after heat treatment will studied at several temperatures (structures of RL and ReL pyCs are similar after 1500°C).

Fig. 1

 1 Fig. 1 Courbes de recuit simules (potentiel total en fonction de la temperature) for several ramps. Simulation times are also noticed (en bas a gauche).

Fig. 2 a

 2 Fig. 2 a) Experimental filtered HRTEM image, b) synthetic HRTEM-like 3D image, c) atomistic model, and d) simulated image of the pyC1 (RL).

Table 1 .

 1 . Structural properties (in percentage) of PyC models constructed from a 16 fringes AP PyC image and a 17 fringes HT PyC image.

	PyC sp	sp 2	sp 3 C 5	C 6	C 7	C 8 pC 6
	AP	0.5 97.3 2.2 5.4 88.6 5.7 0.3 75.3
	HT	0.2 99.0 0.8 5.4 88.9 5.5 0.2 76.4
	Finally, this study showed an excellent agreement
	between experiments and simulations looking at
	High Resolution Transmission Electron Microscopy
	(HRTEM) images and Pair Distribution Functions
	(PDF).

Acknowledgement. Funding from the Agence Nationale de la Recherche through the PyroM aN project (contract ANR-2010-BLAN-929) and from the Institut Carnot M aterials and systems Institute of Bordeaux (M IB) are gratefully acknowledged.


Radu Urs
email: radu.urs@ims-bordeaux.fr

Jean-Pierre Da Costa
email: jean-pierre.dacosta@ims-bordeaux.fr

Jean-Marc Leyssale

Gérard Vignoles

Christian Germain

Gerard Vignoles

Non-parametric synthesis of laminar volumetric texture

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

Many tools in image processing, vision and computer graphics invoke image analysis and synthesis, and especially textured images. For this reason, the disciplinary field of texture synthesis is particularly dynamic and very productive with notable applications in image compression, inpaiting, image extrapolation or texture mapping. 3D textures are mainly used for texturing volumetric objects trying to increase the realism of the 3D scenario. This research field has lead to the development of many 2D synthesis techniques, but their extension to the 3D environment remains unstable proving itself as a very complex and computational issue. However, the interest of the scientific community for solid texture synthesis gave birth to various families of algorithms.

The procedural family of 3D synthesis algorithms models a solid texture by adjusting some mathematical functions. The range of applicability of the procedural approaches [START_REF] Heeger | Pyramid-based texture analysis/synthesis[END_REF][START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF][START_REF] Costa | Synthesis of solid textures based on a 2D example: application to the synthesis of 3D carbon structures observed by transmission electronic microscopy[END_REF] remains sparse, leading mostly to algorithms hard to optimize. In opposition to the global growing character of the solid texture involved by the procedural methods, the local methods [START_REF] Wei | Texture synthesis from multiple sources[END_REF][START_REF] Kopf | Solid Texture Synthesis from 2D Exemplars[END_REF][START_REF] Chen | High quality solid texture synthesis using position and index histogram matching[END_REF][START_REF] Wei | Fast texture synthesis using tree-structured vector quantization[END_REF][START_REF] Paget | Texture synthesis via a noncausal nonparametric multiscale markov random field[END_REF][START_REF] Efros | Texture synthesis by non-parametric sampling[END_REF][START_REF] Lefebvre | Appearance-space texture synthesis[END_REF] generate the texture one voxel/patch at a time maintaining the coherence of the evolving region with its vicinity. Some of the most underlined synthesis methods based on local optimization from the literature are the exemplar-based synthesis algorithms. These algorithms require only one sample texture as input image and synthesize the solid texture guaranteeing the best possible similitude with the sample.

Allowing a smoother control during the synthesis process, the pixel by pixel synthesis algorithms [1-3, 8, 12] are among the most successful techniques mainly because of their practicability. To assure a reasonable computational complexity, most of the approaches propose efficient searching algorithms combined with a multi-scale implementation capable of capturing patterns at different scales and not increasing the computational load.

We are interested in the synthesis of anisotropic volumetric textures with lamellar structure. The synthesis of such laminar textures from a unique sample is however delicate. It means to infer 3D information from 2D information only. This inference is possible only under some assumptions of anisotropy. A recent 2D/3D extension of the Portilla and Simoncelli [START_REF] Portilla | A parametric texture model based on joint statistics of complex wavelet coefficients[END_REF] parametric algorithm was proposed by Da Costa and Germain [START_REF] Costa | Synthesis of solid textures based on a 2D example: application to the synthesis of 3D carbon structures observed by transmission electronic microscopy[END_REF], relying on the 2D to 3D inference of the 1 st and 2 nd order statistics and on an original 3D pyramidal decomposition algorithm. In the present article, we evaluate various implementations of non-parametric algorithms for the synthesis of laminar textures by means of a both qualitative and quantitative procedure. As an application, we address the particular case of carbon material images obtained by HRTEM [START_REF] Leyssale | An image-guided atomistic reconstruction of pyrolytic carbons[END_REF] showing lamellar structure, objective sketched in Fig. 1. Four examples of such laminar textures are provided in Fig. 2.

Non-parametric volumetric synthesis techniques

We are briefly exposing the non-parametric methods we evaluate, namely the algorithms of Wei and Levoy [START_REF] Wei | Texture synthesis from multiple sources[END_REF][START_REF] Wei | Fast texture synthesis using tree-structured vector quantization[END_REF] to serve as a basis, Kopf et al. [START_REF] Kopf | Solid Texture Synthesis from 2D Exemplars[END_REF] -employed for the improvements brought by its non-parametric global texture optimization procedure and finally Chen and Wang [START_REF] Chen | High quality solid texture synthesis using position and index histogram matching[END_REF] -proposing a high quality texture synthesis by integrating the texture optimization with two new kinds of matching histograms.

Overview of Wei and Levoy's algorithm

Based on the Markov field hypothesis, this method relies on texture locality (every pixel is predictable from the few pixels in its neighbourhood) and texture stationarity (spatial statistics invariance by translation), assumptions stated in Wei and Levoy [START_REF] Wei | Fast texture synthesis using tree-structured vector quantization[END_REF]. The 2D/3D extension that we put up here is an adaptation of the algorithm of texture synthesis from multiple sources [START_REF] Wei | Texture synthesis from multiple sources[END_REF]. Starting from a random initialization, we synthesize voxel by voxel by examining the 2D neighbourhood of the current voxel from two orthogonal views of the 3D block (front and side view); this phase implies a TSVQ (Tree-Structured Vector Quantization) search of the best match for each of these 2 neighbourhoods in the same input image. The output pixel corresponds to the average of the two found voxels and we reiterate several times until reaching the same neighbourhoods after 2 consecutive iterations. This 2D/3D extension raises different downsides being constrained by the neighbourhood size and shape or the synthesis scan type. But, by including a multiresolution scheme it proves to be a well adaptable and largely applicable synthesis method.

Overview of Kopf's algorithm

Another issue of Wei and Levoy's algorithm resides in the way of combining the information from the front view with the side view as shown in Fig. 2. The simple averaging technique leads to a loss of dynamics between the original and the synthesized texture. One possible answer is to replace the average with a better combination [START_REF] Kopf | Solid Texture Synthesis from 2D Exemplars[END_REF][START_REF] Kwatra | Graphcut textures: image and video synthesis using graph cuts[END_REF] and, in the same time, adding a histogram matching mechanism [START_REF] Kopf | Solid Texture Synthesis from 2D Exemplars[END_REF][START_REF] Heeger | Pyramid-based texture analysis/synthesis[END_REF]. Kopf et al. [START_REF] Kopf | Solid Texture Synthesis from 2D Exemplars[END_REF] propose to give weights to each 2 found voxels and to mix them according to the robust optimization proposed by Kwatra et al. [START_REF] Kwatra | Graphcut textures: image and video synthesis using graph cuts[END_REF]. The histogram matching works by adjusting the weight of each voxel engaged in the weighted average that could bring differences between the colour histogram of the 3D result and the colour histogram of the 2D exemplar. In this way the global statistics are preserved and the synthesis doesn't rely only on local decisions. The quality of the results delivered with the technique proposed by Kopf exceeds that of many other methods. Despite this, computing the results by using the weighted average may produce blurry results if the variance of the exemplar voxels is too large [START_REF] Kopf | Solid Texture Synthesis from 2D Exemplars[END_REF].

Overview of Chen and Wang's algorithm

Unfortunately, most of the non-parametric methods are affected by blurring, missing textural patterns or mismatched input/output histograms. To ameliorate these drawbacks, Chen and Wang [START_REF] Chen | High quality solid texture synthesis using position and index histogram matching[END_REF] propose a new modus operandi for texture optimization approaches. Mainly, it consists in using two new kinds of histograms -position and index histograms [START_REF] Chen | High quality solid texture synthesis using position and index histogram matching[END_REF] and in updating the output voxel using the discrete solver as Han et al. [START_REF] Han | Fast example-based surface texture synthesis via discrete optimization[END_REF]. Prior to synthesis, for every pixel in the exemplar we have to find the k best matches from the exemplar, idea inspired from the natural texture synthesis algorithm proposed by Ashikmin [START_REF] Ashikhmin | Synthesizing natural textures[END_REF]. The set of candidates being constructed, we start a synthesis process similar to Kopf et al. [START_REF] Kopf | Solid Texture Synthesis from 2D Exemplars[END_REF], but instead of using colour histogram we use the position histogram to adjust the weights. The position histogram records the number of occurrences in the solid texture of the corresponding pixel in the 2D exemplar. The weights are adjusted to get a histogram as uniform as possible. Once the position histogram matching is achieved, the colour histogram matching is achieved also. Index histogram counts the frequency of the voxels candidates in the volume texture. Having all the candidates uniformly distributed means having a well preserved texture. During the search phase, some pixels are restricted according to their index histogram frequency. This algorithm provides results that are similar or in part better than the previous methods.

Whatever the method, a problematic issue remains the choice of the neighbourhood system. The choice of a causal neighbourhood makes the synthesis of a pixel totally dependent on previous pixels and leads in the case of lamellar textures to a higher regularity of synthesized results. The choice of a noncausal neighbourhood can partly overcome this determinism, but now the choice of the scan type may be challenged. Replacing the lexicographical scan type with a completely random walk allows the synthesis of a pixel by freeing itself from its past and so multiplying the possible configurations. However, the convergence time can become prohibitive in this case. To reconcile the deterministic part and the randomness character we propose an alternative scan type, namely the use of space filling curves [START_REF] Sagan | Space-Filling Curves[END_REF] extended to 3D, like Morton code (Z-curve) or fractal curves(e.g. Hilbert curve) illustrated in Fig. 4.

Experimentation background

We were able to implement different variants of the synthesis algorithm and apply them to the volumetric texture synthesis starting from a single 2D exemplar obtained by microscopy.

Fig. 5a and5d show the volumetric textures for Wei and Levoy's approach while Fig. 5b and5e show the results obtained using the backbone of Wei and Levoy but optimized by using the weighted average and the colour histogram matching technique proposed by Kopf et al. [START_REF] Kopf | Solid Texture Synthesis from 2D Exemplars[END_REF]. In Fig. 5c and 5f are represented the 3D blocks obtained by using the discrete solver plus the index and the position histogram matching technique proposed by Chen and Wang [START_REF] Chen | High quality solid texture synthesis using position and index histogram matching[END_REF]. The experimentation framework consists in using 5 pyramidal levels, non-causal neighbourhoods of size 7x7 pixels and, for Chen and Wang, a set of 15 candidates for each pixel. The visual quality of the results is relatively convincing, in terms of dynamics and structure. Measuring the results based on our human perception, we can say that the 3D results and the original images seem to have been created by the same underlying process (see the images from Fig. 2 and Fig. 5 for analogy). The perceptual comparison of the result in Fig. 5d with the ones in Fig. 5e or 5f reveals that, in some cases, the Wei and Levoy's algorithm is not capable to preserve the gray-levels, which is not the case for the algorithms involving histogram adjustments. In addition, taking into account the size of the structures in determining the number of synthesis levels and the neighbourhood size is a critical point in terms of synthesis quality.

Results

Quantitative evaluation

The results interpretation using only our senses remains subjective, so we are constrained to compare them objectively. This original quantitative study evaluates the results taking into account the results gray level dynamics (first order statistics), but also their spatial statistics (autocorrelations) and morphological properties like fringe length, tortuosity and orientation. It aims to evaluate the algorithms capacity to reproduce a 3D texture respecting the statistical properties of the exemplar.

The quantitative indicators from Fig. 6, Fig. 8 and Fig. 9 reveal the 3D textures tendency towards the same structure observed on the HRTEM images. The comparison of the results dynamics from Fig. 6 confirms the visual conclusion attributed to Wei and Levoy's algorithm not being able to conserve contrast because of the averaging operation when combining the two orthogonal views. This objective comparison of the evolution of the 3D blocks statistics towards the ones from the 2D sample strengthens the visual assumptions relative to the significant improvements brought by the histogram matching techniques.

As for the structural properties of the synthesized blocks, we have compared in Fig. 8 the local orientations of the multi-2D block relative to the input sample based on a structure tensor estimated on a 7×7 neighbourhood inspired from Bigun et al. [START_REF] Bigun | Multidimensional orientation estimation with applications to texture analysis and optical flow[END_REF] using image orientation field and its associated orientation histogram as exemplified in Fig. 7. The left plot from Fig. 8 confirms the impression that in some cases the proposed algorithms tend to produce textures more regular than that of the input image, and iterating more doesn't resolve the reproduction of patterns. This phenomenon may be related to the deterministic character of the algorithms, based on local decisions, which tends to produce these very regular or repetitive textures. But this can be also caused by the limited size of the input image which is not large enough to guarantee a sufficient diversity at the output, or the limited size of the output block of which periodicity somehow imposes a certain regularity. The general convergence tendency is also certified by the comparison of the morphological properties. It consists in tracking the texture level curves and describing the elongated patterns contained within the texture. We address the orientation and pattern length measurement techniques as described by Da Costa et al. [START_REF] Costa | Level curve tracking algorithm for textural feature extraction[END_REF] to extract textural features (fringes) and compute their lengths and tortuosity. We count the frequency of appearance of the fringe lengths and classify them in function of their lengths. Fig. 9a shows the frequency of the fringe lengths in the input image in comparison to the frequency of the fringe lengths in the multi-2D synthesized 3D block where fringes were retrieved by manipulating the front slices and the side slices in the block. Similar, we built tortuosity classes and compare the tortuosity values distribution in Fig. 9b.

Conclusions

In this paper, we have brought about several algorithms that achieve with success solid texture synthesis starting from a 2D texture exemplar. Results interpretation was accompanied by an original experimental procedure marking a qualitative and objective evaluation. As an application, we have used these algorithms to the synthesis of anisotropic HRTEM images of dense carbons. Apart from certain limitations, the obtained results are convincing in term of preserving texture gray-levels and texture structure properties.

Figure 1 :

 1 Figure 1: The general schema: obtaining the HRTEM-like snapshots of dense carbons and generating a 3D representation from a single 2D sample.

Figure 2 :

 2 Samples of lattice fringe images of dense pre-graphitic carbons: HRTEM images (a, c) and their filtered versions (b, d). Sample sizes are 128×128 pixels. Resolution is about 5×10 -11 m per pixel.

Figure 3 :

 3 Figure 3: Principle of non-parametric synthesis: extract two neighbourhoods in the output block (front and side view), search the best two neighbourhoods in the input image, combine their central pixel values and modify the output voxel.

Figure 4 :

 4 Illustration of different scanning types, in 2-dimensions (the top column) and their corresponding extensions to 3-dimensions (the bottom column): from left to right, a -lexicographic path, b -random walk, c -Z-curve and d -Hilbert curve.

Figure 5 :

 5 Volumetric results: from top to bottom, the 1 st line represents the 3D views of the solid textures obtained after 10 iterations from the raw HRTEM sample in Fig.2ausing, from left to right, Wei and Levoy (a), Kopf (b), Chen and Wang(c); the 2 nd line (d, e, f) shows the blocks after applying the same algorithms on the filtered sample in Fig. 2d. For c and f, next to the synthetic block are represented, from top to bottom, the images corresponding to the index and position histograms.

Figure 6 :

 6 Plots of different indicators, corresponding to the objective comparison of the 3D textures obtained by synthesizing the HRTEM texture in Fig.2a(the 1 st line) and the sample in Fig.2d(the 2 nd line): gray level histograms (a, d) and 1 st order statistics (b, estandard deviation; c, f -skewness). The procedure consist in comparing input image statistics with the ones obtained after a multi-2D solid blocks analysis (b, c, e, f) or by computing Kullback-Leibler divergence between the histograms of the exemplar and the output block (a, d).

Figure 7 :

 7 Texture orientation: a -the original raw HRTEM sample from Fig. 1a; b -the correspondent orientation field and the colour bar giving the correspondence between colours and orientations; c -the orientation histogram computed from the orientation field.

Figure 8 :

 8 Plots of the spatial structure variation indicator: the Kullback -Leibler divergence between the orientation histograms of the exemplar and the output blocks obtained by applying the algorithms on the exemplar in Fig.1a(the left plot -a) and on the sample in Fig.1d(the right plot -b).

Figure 9 :

 9 Plots of the morphological structure indicators: a -the distribution of fringe lengths in the input raw HRTEM image from Fig. 2a and in the volumetric synthesized textures from Fig. 5a-5c; b -the distribution of the tortuosity values corresponding to the same input image and output blocks used in a.

Acknowledgements

This work has been carried out within the project PyroMaN (Pyrocarbon Matrices at the Nanoscale, http://pyroman.dr15.cnrs.fr) funded by the French National Research Agency (ANR) and by Aerospace Valley, World Competitiveness Cluster in Aeronautics, Space and Embedded Systems of Midi-Pyrenees and Aquitaine. The company Snecma Propulsion Solide (Safran group) has to be also acknowledged for providing the image samples used for the tests.


F Bodenan

F Bourgeois

C Petiot

T Auge

A B Bonfils

C Julcour

F Guyot

J Tremosa

Boukary

A Lassin

P Chiquet

EX-SITU MINERAL CARBONATION: RESOURCES, PROCESS AND ENVIRONMENTAL ASSESSMENTS (CARMEX PROJECT)

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

EX-SITU MINERAL CARBONATION: RESOURCES, PROCESS AND ENVIRONMENTAL ASSESSMENTS (CARMEX PROJECT)

Françoise Bodénan, Florent Bourgeois, Charlotte Petiot, Thierry Augé, Benjamin Bonfils, Carine Julcour-Lebigue, François Guyot, Joachim Tremosa, Aïssa Boukary, Arnault Lassin, et al.

To cite this version:

Françoise Bodénan, Florent Bourgeois, Charlotte Petiot, Thierry Augé, Benjamin Through building a GIS it has been possible to identify potentially favourable areas for developing the sector on a worldwide scale. A first treatment targeted on large amounts of (already crushed) ultramafic rock waste near thermal power plants has enabled us to determine the bottle-necks, notably concerning the volumes that are actually accessible, and to identify relevant selection criteria (1). The case of South Africa and, to a lesser extent, that of New Caledonia in connection with its insularity, are discussed. Indeed, New Caledonia's insularity and its accessible abundance of suitable rocks and industrial waste (i.e. rich in MgO, CaO, if not Fe(II)O) would appear to be perfect for the application of mineral carbonation to store the extensive industrial emissions of CO 2 derived from pyrometallurgical processing of nickel.

A detailed study of direct aqueous mineral carbonation in the presence of organic ligands (3) was studied in detail as part of a thesis (2), but detailed analysis of the coupled dissolution and precipitation mechanisms has led this option being discarded (4). Very promising yields (70-90% conversion) under limited P and T conditions, on the other hand, have been obtained by integrating an appropriate physical process.

Comparative environmental assessments of three scenarios for a coal-fired power stationwithout CO 2 capture, with CO 2 capture and geological storage, and with CO 2 capture and exsitu mineral carbonationshow mixed results with pollution gains or transfers depending on the process and the studied impact indicators. The analysis shows that the viability of this CO 2 storage option is mainly at the level of the mineral carbonation process and the optimization of its operating conditions.

This research (Carmex Project No. ANR-08-PCO2-002) is co-funded by the French National Research Agency under the CO 2 Capture and Storage programme and by Total (thesis). Special thanks to P. Galle-Cavalloni for the experiments and to the BRGM laboratory services.