Romain Deville

Elisa Fromont

Baptiste Jeudy

Christine Solnon

GriMa: a Grid Mining Algorithm for Bag-of-Grid-Based Classification

General-purpose exhaustive graph mining algorithms have seldom been used in real life contexts due to the high complexity of the process that is mostly based on costly isomorphism tests and countless expansion possibilities. In this paper, we explain how to exploit gridbased representations of problems to efficiently extract frequent grid subgraphs and create Bag-of-Grids which can be used as new features for classification purposes. We provide an efficient grid mining algorithm called GriMA which is designed to scale to large amount of data. We apply our algorithm on image classification problems where typical Bagof-Visual-Words-based techniques are used. However, those techniques make use of limited spatial information in the image which could be beneficial to obtain more discriminative features. Experiments on different datasets show that our algorithm is efficient and that adding the structure may greatly help the image classification process.

Introduction

General-purpose exhaustive graph mining algorithms are seldom used in realworld applications due to the high complexity of the mining process mostly based on isomorphism tests and countless expansion possibilities during the search [START_REF] Cook | Mining Graph Data[END_REF]. In this paper we tackle the problem of exhaustive graph mining for grid graphs which are graphs with fixed topologies and such that each vertex degree is constant, e.g., 4 in a 2D square grid. This grid structure is naturally present in many boardgames (Checkers, Chess, Go, etc.) or to model ecosystems using cellular automata [START_REF] Hogeweg | Cellular automata as a paradigm for ecological modelling[END_REF], for example. In addition, this grid structure may be useful to capture low-level topological relationships whenever a high-level graph structure is not obvious to design. In computer vision in particular, it is now widely acknowledged that high-level graph-based image representations (such as region adjacency graphs or interest point triangulations, for example) are sensitive to noise so that slightly different images may result in very different graphs [START_REF] Prado | Mining spatiotemporal patterns in dynamic plane graphs[END_REF][START_REF] Ozdemir | Image classification using subgraph histogram representation[END_REF]. However, at a low-level, images basically are grids: 2D grids of pixels for images, and 3D grids of voxels for tomographic or MRI images modelling 3D objects. When considering videos, we may add a temporal dimension to obtain 2D+t grids. We propose to exploit this grid structure and to characterize images by histograms of frequent subgrid patterns.

Section 2 describes existing works that use pattern mining approaches for image classification and existing graph mining algorithms related to our proposed approach. Section 3 introduces grid graphs and our sub-grid mining algorithm GriMA. Section 4 experimentally compares these algorithms and shows the relevance of GriMA for image classification.

Related Work

Mining for image classification. Pattern mining techniques have recently been very successfully used in image classification [START_REF] Fernando | Mining mid-level features for image classification[END_REF][START_REF] Silva | Bog: A new approach for graph matching[END_REF] as a mean to obtain more discriminative mid-level features. However, those approaches consider the extracted features used to describe images (e.g., Bags-of-Visual-Words / BoWs) as spatially independent from each other. The problem of using bag-of-graphs instead of BoWs has already been mentioned in [START_REF] Da Silva | Image classification based on bag of visual graphs[END_REF][START_REF] Odone | Building kernels from binary strings for image matching[END_REF][START_REF] Samuel | Extracting plane graphs from images[END_REF] for satellite image classification and biological applications. However, none of these papers provide a general graph representation nor a graph mining algorithm to extract the patterns. In [START_REF] Nowozin | Weighted substructure mining for image analysis[END_REF], authors have already shown that, by combining graph mining and boosting, they can obtain classification rules based on subgraph features that contain more information than sets of features. The gSpan algorithm [START_REF] Voravuthikunchai | Histograms of Pattern Sets for Image Classification and Object Recognition[END_REF] is then used to compute the subgraph patterns but a limited number of features per image is used to be able to scale on real-life datasets.

Graph Mining. gSpan [START_REF] Voravuthikunchai | Histograms of Pattern Sets for Image Classification and Object Recognition[END_REF] and all similar general exhaustive graph mining algorithms [START_REF] Jiang | A survey of frequent subgraph mining algorithms[END_REF] extract frequent subgraphs from a base of graphs. During the mining process, gSpan does not consider edge angles so that it considers as isomorphic two subgraphs that are different in a grid point of view as shown in Fig. 1 (b andc). Because of this, gSpan may consider as frequent a huge number of patterns and does not scale well. On the other hand, Plagram [START_REF] Ozdemir | Image classification using subgraph histogram representation[END_REF] has been developed to mine plane graphs and thus to scale efficiently on these graphs. However, in Plagram, the extension strategy (which is a necessary step for all exhaustive graph mining algorithms) is based on the faces of the graph which induces an important restriction: All patterns should be composed of faces and the smallest possible subgraph pattern is a single face, i.e, a cycle with 3 nodes. Using Plagram to mine grids is possible but the problem needs to be transformed such that each grid node becomes a face. This transformation is illustrated in Fig. 1: The two isomorphic graphs (b) and (c) (which are both subgraphs of (a)) are transformed into the two non-isomorphic graphs (e) and (f) (so that only (e) is a subgraph of (d)). However, this artificially increases the number of nodes and edges which may again cause scalability issues. Grid graphs have already been introduced in [START_REF] Marinescu-Ghemeci | Maximum induced matchings in grids[END_REF], but authors did not take into account the rigid topology of the grid and, in particular, the angle information that we use to define our grid mining problem. Finally, grid graphs are special cases of geometric graphs, for which mining algorithms have been proposed in [START_REF] Arimura | Time and space efficient discovery of maximal geometric graphs[END_REF]. These algorithms (FreqGeo and MaxGeo) may be seen as a generalization of our new algorithm GriMA but, as such, they are not optimized for cases where the graph is known to be a grid and we show in Sec. 3 that their complexity is higher than the complexity of our algorithm. Besides, the authors do not provide any implementation (thus, also no experiment) of their proposed algorithm which could allow us a comparison with our method.

Grid Mining

Grid Graphs. A grid is a particular case of graph such that edges only connect nodes which are neighbors on a grid. More formally, a grid is defined by G = (N, E, L, η) such that N is a set of nodes, E ⊆ N × N is a set of edges, L : N ∪E → N is a labelling function, η : N → Z 2 maps nodes to 2D coordinates, and ∀(u, v) ∈ E, the coordinates of u and v are neighbors, i.e., |x u -x v |+|y u -y v | = 1 where η(u) = (x u , y u ) and η(v) = (x v , y v ). Looking for patterns in a grid amounts to searching for subgrid isomorphisms. In an image analysis context, patterns should be invariant to translations and rotations so that grid patterns which are equivalent up to a translation and/or a rotation should be isomorphic. For example, graphs (b) and (c) of Fig. 1 are isomorphic. However, (c) cannot be obtained by translating and/or rotating (b), because the angle between edges of (c) is different from the angle between edges of (b). Therefore, (b) and (c) are not grid-isomorphic. Finally, a grid

G 1 = (N 1 , E 1 , L 1 , η 1 ) is subgrid-isomorphic to a grid G 2 = (N 2 , E 2 , L 2 , η 2 ) if there exists N 2 ⊆ N 2 and E 2 ⊆ E 2 ∩ N 2 × N 2 such that G 1 is grid isomorphic to G 2 = (N 2 , E 2 , L 2 , η 2 ).
Frequent Grid Mining Problem. Given a database D of grid graphs and a frequency threshold σ, the goal is to output all frequent subgrid patterns in D, i.e., all grids G such that there exist at least σ grids in D to which G is subgrid isomorphic.

Grid Codes. In graph mining algorithms, an important problem is to avoid generating the same pattern multiple times. One successful way to cope with this is to use canonical codes [START_REF] Voravuthikunchai | Histograms of Pattern Sets for Image Classification and Object Recognition[END_REF] to represent graphs and thus explore a canonical code search space instead of a graph one. In this paragraph, we define the canonical code used in our algorithm. A code C(G) of a grid G is a sequence of n edge codes (C(G) = ec 0 , ..., ec n-1 ) which is associated with a depth-first traversal of G starting from a given initial node. During this traversal, each edge is traversed once, and nodes are numbered: The initial node has number 0; Each time a new node is discovered, it is numbered with the smallest integer not already used in the traversal. Each edge code corresponds to a different edge of G and the order of edge codes in C(G) corresponds to the order edges are traversed. Hence, ec k is the code associated with the k th traversed edge. This edge code ec k is the tuple (δ, i, j, a, L i , L j , L (i,j) ) where i and j are the numbers associated with the nodes of the k th traversed edge.

δ ∈ {0, 1} is the direction of the k th traversed edge:

• δ = 0 if it is forward, i.e., i already appears in the prefix ec 0 , ..., ec k-1 of the code and j is a new node which is reached for the first time; • δ = 1 if it is backward, i.e., both i and j already appear in the prefix.

a ∈ {0, 1, 2, 3} is the angle of the k th traversed edge:

• a = 0 if k = 0 (first edge);

• Otherwise, a = 2A/π where A ∈ {π/2, π, 3π/2} is the angle between the edge which has been used to reach i for the first time and (i, j). -L i , L j , L (i,j) are node and edge labels.

For example, let us consider code 1 in Fig. 2. The fourth traversed edge is (c, d). It is a forward edge (because d has not been reached before) so that δ = 0. The angle between (b, c) and (c, d) is 3π/2 so that a = 3. The fifth traversed edge is (d, a) which is a backward edge (as a has been reached before) so that δ = 1.

Given a code, we can reconstruct the corresponding grid since edges are listed in the code together with angles and labels. However, there exist different possible codes for a given grid, as illustrated in Fig. 2: Each code corresponds to a different traversal (starting from a different initial node and choosing edges in a different order). We define a total order on the set of all possible codes that may be associated with a given grid by considering a lexicographic order (all edge code components have integer values). Among all the possible codes for a grid, the largest one according to this order is called the canonical code of this grid and it is unique. Description of GriMA. Our algorithm, GriMA, follows a standard frequent subgraph mining procedure as described, e.g., in [START_REF] Ozdemir | Image classification using subgraph histogram representation[END_REF]. It explores the search space of all canonical codes in a depth-first recursive way. The algorithm first computes all frequent edges and then calls an Extend function for each of these frequent extensions. Extend has one input parameter: A pattern code P which is frequent and canonical. It outputs all frequent canonical codes P such that P is a prefix of P . To this aim, it first computes the set E of all possible valid extensions of all occurrences of P in the database D of grids: A valid extension is the code e of an edge such that P.e occurs in D. Finally, Extend is recursively called for each extension e such that P.e is frequent and canonical. Hence, at each recursive call, the pattern grows.

Node-induced GriMA. In our application, nodes are labelled but not edges (all edges have the same label). Thus, we designed a variant of GriMA, called node-induced-GriMA, which computes node-induced grids, i.e., grids induced by their node sets. This corresponds to a "node-induced" closure operator on graphs where, given a pattern P , we add all possible edges to P without adding new nodes. In the experiments, we show that this optimization decreases the number of extracted patterns and the extraction time.

Properties of GriMA and Complexity. We can prove (not detailed here for lack of space) that GriMA is both correct, which means that it only outputs frequent subgrids and complete, which means that it cannot miss any frequent subgrid. Let k be the number of grids in the set D of input grids, n the number of edges in the largest grid G i ∈ D and |P | the number of edges in a pattern P . GriMA enumerates all frequent patterns in O(kn 2 .|P | 2 ) = O(kn 4 ) time per pattern P . This is a significant improvement over FreqGeo and MaxGeo [START_REF] Arimura | Time and space efficient discovery of maximal geometric graphs[END_REF] which have a time complexity of O(k 2 n 4 . ln n) per pattern.

Experiments

To assess the relevance of using a grid structure for image classification, we propose to model images by means of grids of visual words and to extract frequent subgrids to obtain Bags-of-Grids (BoGs). We compare these BoGs with a standard classification method which uses simple unstructured Bags-of-Visual-Words (BoWs). Note that neither BoG nor BoW-based image classification give stateof-the-art results for these datasets. In particular, [START_REF] Chatfield | Return of the devil in the details: Delving deep into convolutional nets[END_REF] reported that the features discovered using deep learning techniques give much better accuracy results than BoWs and all their extensions on classification problems. However, our aim is to compare an unstructured set of descriptors and a set of descriptors structured by the grid topology. The method presented in this paper is generic and may be used with any low-level features (e.g. deep-learned features) as labels.

Datasets. We consider three datasets. Flowers [START_REF] Nilsback | Automated flower classification over a large number of classes[END_REF] contains 17 flower classes where each class contains 80 images. 15-Scenes [START_REF] Lazebnik | Beyond bags of features: Spatial pyramid matching for recognizing natural scene categories[END_REF] contains 4485 images from 15 scene classes (e.g., kitchen, street, forest, etc.) with 210 to 410 images per class. Caltech-101 [START_REF] Fei-Fei | Learning generative visual models from few training examples: An incremental bayesian approach tested on 101 object categories[END_REF] contains pictures of objects belonging to 101 classes. There are 40 to 800 images per class. As it is meaningless to mine frequent patterns in small datasets of images (when there are not enough images in a dataset, there may be a huge number of frequent patterns), we consider a subset of Caltech-101 which is composed of the 26 classes of Caltech-101 that contain at least 80 pictures per class. For each class, we randomly select 80 images. This subset will be referred to as Caltech-26 BoW design. Given a set of patches (small image regions) extracted from an image dataset, visual words are created by quantizing the values of patch descriptors [START_REF] Csurka | Visual categorization with bags of keypoints[END_REF][START_REF] Chatfield | The devil is in the details: an evaluation of recent feature encoding methods[END_REF]. The set of computed visual words is called the visual vocabulary. Each image is then encoded as an histogram of this visual vocabulary, called a bag-of-visual-words. A decade ago, patches used to create the visual vocabulary were selected by using interest point detectors or segmentation methods. However, [START_REF] Nowak | Sampling strategies for bag-of-features image classification[END_REF] has shown that randomly sampling patches on grids (called dense sampling) gave as good (and often better) results for image classification than when using complex detectors. This made this problem a particularly suited use-case for our algorithm. In our experiments, visual words are based on 16x16 SIFT descriptors [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF] which are 128-D vectors describing gradient information in a patch centered around a given pixel. 16x16 SIFT descriptors are extracted regularly on a grid and the center of each SIFT descriptor is separated by s = 8 pixels (descriptors are thus overlapping). We use the K-means algorithm to create visual word dictionaries for each dataset. The optimal number of words K is different from one method (BoG) to the other (BoW) so this parameter is studied in the experiments.

BoG design. The first steps (computation of descriptors densely distributed in each image and creation of the visual vocabulary) are similar for BoW and BoGbased methods. For BoG, we create a square grid based on the grid of visual words by connecting each node to its 4 neighbors (except on the border of the image). In our experiments, grid nodes are labeled by visual words and edges remain unlabeled although the algorithm is generic and could include labels on edges. For efficiency reasons, we preprocess grids before the mining step to remove all nodes which are surrounded by nodes with identical labels (thus creating holes in place of regions with uniform labels). Once grids are created for all images, we mine frequent subgrid patterns with GriMA class by class. Finally, we consider the union of all frequent patterns and all visual words, and represent each image by an histogram of frequency of patterns and words. These histograms are given as input to the classifier. Efficiency Results. Figure 3 compares scale-up properties of GriMA, nodeinduced-GriMA, gSpan and Plagram for the Flowers dataset (average results on 10 folds for each class, where each fold contains 72 images). We have set the number K of words to 100, and we have increased the absolute frequency threshold from 45 (62.5%) to 72 (100%). Note that, as explained in Sect. 1. Classification results according to the number of visual words K and the datasets. For BOG, the column #pat displays the number of frequent patterns extracted by node-induced-GriMA and σ give the frequency freshold (in percentage of the number of images per class) reached to obtain at least 4000 patterns.

the size K of the vocabulary and on the dataset. Therefore, instead of fixing σ, we propose to fix the number p of patterns we want to mine: Starting with σ = 100%, we iteratively decrease σ and run node-induced-GriMA for each value of σ until the number of mined patterns reaches p. In our experiments, we have set p to 4, 000 as this value corresponds to a good tradeoff between time and gain in accuracy.

For the classification task, we use N (N -1)/2 binary Support Vector Machines (SVM), where N is the number of classes. Each SVM is trained to discriminate one class against another. The final classification is obtained using a vote among the classes. SVMs are trained using Libsvm [START_REF] Chang | LIBSVM: A library for support vector machines[END_REF] with the intersection kernel presented in [?]. Inputs of SVMs are either BoW histograms (for the unstructured visual word-based method) or BoW ∪ BoG histograms (for the grid-based method). The parameter C of SVM has been optimized by cross validation on all datasets. We consider 10-fold cross-validation for Flowers and Caltech-26 since the number of images is the same for each class. For 15-Scenes, we have created 10 different folds: For each fold, we randomly select 100 images per class for training and 50 images per class for testing. We provide a Student two-tailed t-test to statistically compare both accuracies (BoW and BoW ∪ BoG) for each dataset (the results in green are statistically better, in red statistically worse at p = 0.05 level).

We report in Table 1 accuracy results obtained for all datasets. As expected, we need to lower the support threshold σ to obtain at least 4000 patterns when the number of visual words increases (since there are more different labels in the grid, the patterns are less frequent). On 15-Scenes, the use of structural information becomes important only for K = 4000 (but this also gives the best accuracy results). On Flowers on the contrary, BOG greatly improves the results for low number of visual words ([100,1000]) but it is not statistically better (nor worse) for higher numbers ([2000,4000]) where the accuracy is the best (K = 2000). For Caltech-26, BOG slightly improves the results when K = 100, and it is not significantly different for higher values of K. Overall, those results show that adding structural information does not harm and sometimes improve the representation of images and that our grid mining algorithm can be used in a real-life context. Binary classification results. Finally, we evaluate the interest of using grids for a binary classification task on the 15-scenes dataset (similar results were observed for the two other datasets). The goal is to provide an insight into each class separately, to see if all classes benefit from the use of grids in a same way. Table 2 reports binary classification results: For each class C, we train a binary SVM to recognize images of C, against images of all other classes. We consider the same parameter settings as for the N -ary classification task described previously, except that we train N binary SVMs (instead of N (N -1)/2). Only the 7 classes with statistically significant differences are shown in the table. We can see that some classes really benefit from the use of structured patterns for almost all frequency tresholds (e.g., industrial, coast, suburb, forest) whereas for some classes, using unstructured information gives better results (e.g., street or moutain). This is due to the fact that for some classes, the structure is too similar from this class to another to use it to discriminate classes.

Conclusion

We have presented GriMA, a grid mining algorithm for Bag-of-Grid-based classification. GriMA takes into account relative angle information between nodes in the grid graph (as well as possible labels on the nodes and the edges) and is designed to scale to large amount of data. We believe that the grid structure can not only be useful for directly modeling some real life problems but can also be used when an underlying high level structure is not obvious to define as it is the case for image classification problems. Experiments show that our algorithm is necessary (compared to existing ones) to efficiently tackle real life problems. They also show on three image datasets that patterns extracted on this grid structure can improve the classification accuracies. However, to further increase the discriminative power of the grid-patterns for image classification we would need to combine at the same time state-of-the-art deep-learned descriptors and design smart post-processing steps as the ones developed in [START_REF] Fernando | Mining mid-level features for image classification[END_REF] for unstructured models. Besides, we plan to upgrade the GriMA algorithm to mine 2D+t grids which is necessary to tackle different real-life applications such as the analysis of ecosystems.

Fig. 1 .

 1 Fig. 1. Examples. (b) and (c) are isomorphic graphs and are both subgraphs of (a). However, (b) and (c) are not grid isomorphic, and (b) is a subgrid of (a) whereas (c) is not a subgrid of (a). (d) (resp. (e) and (f)) is obtained from (a) (resp. (b) and (c)) by replacing each node by a 4-node face (with the same label on the 4 nodes). (e) and (f) are not isomorphic, and (e) is a subgraph of (d) whereas (f) is not a subgraph of (d).

Fig. 2 .

 2 Fig. 2. Examples of grid codes (node labels are displayed in blue; edges all have the same label 0). Codes 1 and 2 correspond to traversals started from edge (a, b), and differ on the third edge. Code 3 corresponds to a traversal started from edge (d, c) and it is canonical.

Table 2 .

 2 Accuracy of each binary SVM for 7 classes of the 15-scenes dataset with a vocabulary of K = 1000 visual words and a threshold σ ∈ [75%, 20%]. For each class, we highlight in bold BOG results that are better than BOW, and in green (resp. red) BOG results that are significantly better (resp. worse) than BOW. For the other 8 classes, BOW and BOG are never significantly different.

		BOW	BOG 75% 70% 65% 60% 55% 50% 45% 40% 35% 30% 25% 20%
	suburb	99.1 99.1 99.1 99.1 99.1 99.2 99.2 99.2 99.2 99.3 99.3 99.3 99.3
	coast	97.4 97.4 97.5 97.5 97.5 97.5 97.5 97.5 97.5 97.5 97.2 97.1 97.2
	forest	98.8 98.9 98.9 98.9 98.9 98.9 98.9 98.9 98.9 98.9 98.9 98.8 98.7
	highway	97.5 97.7 97.7 97.6 97.7 97.6 97.5 97.5 97.4 97.3 97.3 97.2 97.1
	moutain	97.7 97.6 97.6 97.6 97.7 97.6 97.6 97.5 97.5 97.5 97.5 97.4 97.3
	street	97.8 97.9 97.8 97.8 97.7 97.7 97.8 97.7 97.7 97.5 97.3 97.1 96.8
	industrial 94.6 94.7 94.9 94.8 95.0 95.1 95.1 95.1 95.2 95.3 95.6 95.4 95.4

Acknowledgement

This work has been supported by the ANR project SoLStiCe (ANR-13-BS02-0002-01).

Code 1

Code 2 Code 3 edge δ i j a L i L j L ij edge δ i j a L i L j L ij edge δ i j a L i L j L ij (a, b) 0 0 1 0 0 1 0 (a, b) 0 0 1 0 0 1 0 (d, c) 0 0 1 0 2 1 0 (b, c) 0 1 2 3 1 1 0 (b, c) 0 1 2 3 1 1 0 (c, e) 0 1 2 2 1 1 0 (c, e) 0 2 [START_REF] Chang | LIBSVM: A library for support vector machines[END_REF] 1 1 1 0 (c, d) 0 2 3 3 1 2 0 (c, b) 0 1 3 1 1 1 0 (c, d) 0 2 4 3 1 2 0 (d, a) 1 3 0 3 2 0 0 (b, a) 0 3 4 1 1 0 0 (d, a) 1 4 0 3 2 0 0 (c, e) 0 2 4 1 1 1 0 (a, d) 1 4 0 1 0 2 0

face-based expansion strategy of Plagram does not allow it to find patterns with no face (such as trees, for example). To allow Plagram to mine the same patterns as its competitors, each node of the original grid graph is replaced by a four-node face (as illustrated in Fig. 1). This way, each frequent subgrid in the original grid is found by Plagram in the expanded grid. However, some patterns found by Plagram do not correspond to patterns in the original grid (e.g., faces in the expanded grid which correspond to edges without their end nodes in the original grid). For this reason, the number of patterns that are mined by Plagram and its computation time are higher than the ones reported for its competitors. gSpan does not consider edge angles when mining subgraphs, and two isomorphic subgraphs may have different edge angles so that they do not correspond to isomorphic subgrids, as illustrated in Fig. 1. Therefore, gSpan and GriMA compute different sets of frequent patterns. However, Fig. 3 shows us that the number of patterns is rather similar (slightly higher for gSpan than for GriMA). GriMA and node-induced-GriMA scale better than gSpan, and are able to find all frequent patterns in less than 10s when σ = 60 whereas gSpan has not completed its execution after 3600s. The number of patterns found by both Plagram and gSpan for a support lower than 65 are thus not relevant and are not shown in the Figure . As explained in Sect. 3, node-induced-GriMA mines less patterns and is more efficient than GriMA, and much more efficient than gSpan and Plagram.

N -ary classification results. Let us now evaluate the interest of using grids for image classification on our three datasets. For each dataset we consider different values for the number K of visual words, with K ∈ {100, 500, 1000, 2000, 4000}. Note that this corresponds to changing the K of K-means.

GriMA has one parameter, σ, and Fig. 3 shows us that the number of mined patterns grows exponentially when decreasing σ. CPU time is closely related to the number of mined patterns, and σ should be set so that we have a good tradeoff between the time needed to mine patterns and the gain in accuracy. The number of mined patterns (and the CPU time) depends on σ, but also on


Gabriel Rilling
email: grilling@ens-lyon.fr.

Fellow, IEEE Patrick Flandrin
email: flandrin@ens-lyon.fr.

Paulo Gonçalves
email: paulo.goncalves@inria.fr.

Jonathan M Lilly
email: jonathanlilly@gmail.com

Bivariate Empirical Mode Decomposition

Keywords: Empirical Mode Decomposition, complex-valued signals, bivariate time series EDICS Category: DSP-TFSR I

. The method being initially limited to real-valued time series, we propose here an extension to bivariate (or complex-valued) time series which generalizes the rationale underlying the EMD to the bivariate framework. Where the EMD extracts zero-mean oscillating components, the proposed bivariate extension is designed to extract zero-mean rotating components. The method is illustrated on a real-world signal and properties of the output components are discussed.

coordinates of a point moving in a 2-dimensional space. In particular, the meaning of the signal should not depend on the choice of such Cartesian coodinates. It is worth noticing that another bivariate extension has been introduced very recently [START_REF] Tanaka | Complex Empirical Mode Decomposition[END_REF]. The difference with the one we propose here is significant, but in a nutshell that other method cleverly uses the original EMD to decompose bivariate time series, while ours is a new algorithm that adapts the rationale underlying the EMD to the bivariate framework. Further comparison of the approaches however is out of the scope of this paper. The communication is organized as follows. The bivariate extension is introduced in Section II. Section III is about the components of the resulting decomposition and an illustration is proposed in Section IV. Additionally, free Matlab/C codes corresponding to the proposed algorithms are made available at http://perso.ens-lyon. fr/patrick.flandrin along with small scripts aimed at reproducing the figures and other EMDrelated software.

II. FROM UNIVARIATE EMD TO BIVARIATE EMD

A. Classical EMD

Basically, the EMD considers a signal at the scale of its local oscillations. The main idea of EMD is then to formalize the idea that, locally: "signal = fast oscillations superimposed on slow oscillations".

Looking at a single oscillation (defined, e.g., as the signal between two consecutive local minima), the EMD is designed to define a local "low frequency" component as the local trend m 1 [x](t), supporting a local "high frequency" component as a zero-mean oscillation or local detail d 1 [x](t), so that we can express x(t) as

x(t) = m 1 [x](t) + d 1 [x](t). (1)
By construction, d 1 [x](t) is an oscillatory signal and, if it is furthermore required to be locally zero-mean everywhere, it corresponds to what is referred to as an Intrinsic Mode Function (IMF) [START_REF] Huang | The Empirical Mode Decomposition and Hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF]. Practically, this primarily implies that all its maxima are positive and all its minima are negative. On the other hand, all we know about m 1 [x](t) is that it locally oscillates more slowly than d 1 [x](t). We can then apply the same decomposition to it, leading to

m 1 [x](t) = m 2 [x](t) + d 2 [x](t)
and, recursively applying this on the m k [x](t), we get a representation of x(t) of the form

x(t) = m K [x](t) + K k=1 d k [x](t). (2)
The discrimination between "fast" and "slow" oscillations is obtained through an algorithm referred to as the sifting process [START_REF] Huang | The Empirical Mode Decomposition and Hilbert spectrum for nonlinear and non-stationary time series analysis[END_REF] which iterates a nonlinear elementary operator S on the signal until some stopping criterion is met. Given a signal x(t), the operator S is defined by the following procedure: Identify all extrema of x(t)

1
Interpolate (using a cubic spline) between minima (resp.

2 maxima), ending up with some "envelope" e min (t) (resp. e max (t))

Compute the mean m(t) = (emin(t)+emax(t))

2 3 Subtract [to] from the signal to obtain S[x](t) = x(t) -m(t) 4
If the convergence criterion is met after n iterations, the local detail and the local trend are defined

as d 1 [x](t) = S n [x](t) and m 1 [x](t) = x(t) -d 1 [x](t).

B. Envelopes in 3 dimensions

The EMD is based on the intuitive notion of "oscillation" which naturally relates to local extrema. But the notion of oscillation is much more confusing when the analyzed data is intrinsically bivariate and it is unclear how to define and interpret local extrema. What is rather clear on the other hand is the notion of rotation, which moreover is arguably a two-dimensional extension of the usual notion of a univariate oscillation. Therefore, the basic idea underlying the proposed bivariate EMD is to formalize the following idea: "bivariate signal = fast rotations superimposed on slower rotations". As with the classical EMD, it is clear that the adopted viewpoint is a priori rather restrictive as, e.g., a white noise signal is not meaningfully treated as a sum of oscillations (or rotations). Still this does not prevent the algorithm from producing a decomposition for any signal, as with the univariate EMD.

In order to separate the more rapidly rotating component from slower ones, the idea is once again to define the slowly rotating component as the mean of some "envelope". Yet the envelope is now a 3-dimensional tube that tightly encloses the signal (see Fig. 1 (b)). Given this, the slowly rotating portion of the signal at any point in time can then be defined as the center of the enclosing tube. To this end, only a given number of points on the tube's periphery are considered, each one being associated with a specific direction. If only 4 points are used, these can be the extreme points in the directions top, bottom, left and right (see Fig. 2). In practice, the top point, for example, is uniquely defined only when the signal reaches a local maximum in the vertical direction and is therefore tangent to the top of the tube. Between such characteristic moments in time, the top point is then simply defined using interpolation, ending up with the "envelope" associated with the upwards direction (cf the black thick lines in Fig. 1 (b)). Now, given some set of points on the tube periphery at a given instant in time, there are at least two ways to define their mean:

1) define the mean as the barycenter of the 4 points, considering each to have unit mass (see Fig. 2 (a)).

2) define the mean as the intersection of two straight lines, one being halfway between the two horizontal tangents, the other one halfway between the vertical ones (see Fig. 2 (b)).

In practice, however, the second scheme may be preferred because it is naturally more robust to sampling errors. More precisely, the reason for this is that the envelope points are defined up to an uncertainty that is not isotropic. Indeed, the order of magnitude of this uncertainty can be estimated through a Taylor expansion, which results in an uncertainty jointly proportional to dx/dt and to the sampling period. Thus, the uncertainty is greatest in the direction locally tangent to the signal, and much smaller (of second order) in the orthogonal direction. As the second scheme only uses information from the orthogonal direction, it is naturally more accurate, especially when the signal is sampled sparsely with respect to its period. Note that sampling effects shall not be taken too lightly as the original EMD has been shown to be very sensitive to sampling [START_REF] Rilling | On the influence of sampling on the Empirical Mode Decomposition[END_REF].

The desired goal concerning the interpolation is the same as for the classical EMD: a smooth interpolation with as few "spurious bumps" as possible. Among common interpolation schemes, this calls for cubic spline as it is well known for its minimum curvature property and, in practice, it is still considered the best interpolation scheme for the EMD [START_REF] Huang | Hilbert-Huang Transform and Its Applications[END_REF].

In the preceding discussion, we have limited ourselves to 4 directions for the sake of simplicity, but there is of course no such limitation in practice. Moreover a large number of directions may be interesting insofar as it reduces the dependance of the final decomposition with respect to rotations of the spatial coordinates. For convenience of the presentation, the bivariate time series in the following will be treated as complex-valued time series. Given a set of directions ϕ k = 2kπ/N, 1 ≤ k ≤ N , the proposed bivariate extensions are defined by the same algorithm as the basic EMD, only with new sifting elementary operators S B1 and S B2 corresponding respectively to the algorithms Algo. 1 and Algo. 2.

Algorithm 1: EMD bivariate extension: scheme 1

for 1 ≤ k ≤ N do 1
Project the complex-valued signal x(t) on direction ϕ k :

2 p ϕk (t) = Re e -iϕk x(t)
Extract the locations {t k i } of the maxima of p ϕk (t).

3

Interpolate the set {(t k i , x(t k i ))} to obtain the 4 envelope curve in direction ϕ k : e ϕk (t).

Compute the mean of all envelope curves:

m(t) = 1 N k e ϕk (t) 5
Subtract the mean to obtain

S B1 [x](t) = x(t) -m(t) 6
Furthermore, the second algorithm can be greatly simplified if we notice that the set that is interpolated at step 4 is in fact included in the plane containing the time axis and direction ϕ k . Thus, the interpolation that is performed is very similar to that in the original EMD sifting. Hence, if the number of considered

S B1 x(t) ≈ x(t) or S B2 x(t) ≈ x(t). (3)

Time (days)

Fig. 3. A signal and its Bivariate Empirical Mode Decomposition. The real parts are plotted as solid blue lines and the imaginary parts as dashed black lines. The decomposition has been obtained with the second algorithm using 64 directions and 10 iterations to extract each component. Portions where components are rotating can be identified by a constant phase shift between the real and imaginary parts: if these are in quadrature, then the signal is rotating circularly; other phase shifts correspond to elliptic rotations.

For the sake of simplicity, we will only address the case of simple periodic solutions which are convenient to describe while still being rather general. Indeed, the method operating at a local scale, its behavior on signals whose properties evolve slowly with respect to the local period is very similar to its operation on exactly periodic signals of constant amplitude. Moreover, signals will here also be considered infinitely continuously differentiable, since the purpose of this section is to provide insight into some possible solutions rather than to make an exhaustive study.

What we find from various simulations is that both algorithms generally accept two types of solutions, the first one corresponding to rotating signals as intended and the second one to cases where the method actually fails to extract rotating components and therefore outputs signals that wander around zero in a more complicated way. It is worth noting that examples of the second type of solutions are generally encountered when the analyzed signal does not clearly contain rotating components, as in e.g. a complexvalued white gaussian noise signal. Moreover, it seems that these non-rotating solutions almost never occur when the analyzed signal is analytic (or anti-analytic). On the other hand, solutions from the first type are signals in which the local sense of rotation never changes. The latter can be defined, e.g., from the sense of the vector product of the velocity and acceleration vectors:

sign Im dx dt • d 2 x dt 2 * . (4)
Alternatively, if the signal rotates in a conterclockwise fashion (and |dx/dt| > 0, ∀t), this exactly means that the signal's derivative can be expressed as dx dt (t) = r(t)e iψ(t) , r(t) > 0, dψ(t) dt > 0.

(

Finally, periodic signals satisfying [START_REF] Richardson | Tracking three Meddies with SOFAR floats[END_REF] are not necessarily fixed points of the sifting operators (3). One of the reasons for this is that there is no condition specifying that they are rotating around zero yet.

In order to clarify what is meant by "around zero", we can consider the very simple case in which the signal performs only one full rotation around zero per period (ψ(t + T ) = ψ(t) + 2π where T is the period). This property simplifies greatly the study, as one can easily show that it implies there is only one maximum per period at step 2 in the algorithms. Therefore, all envelope curves are in fact constants with respect to time, which allows one to derive the mean analytically. The envelope curve associated with the direction ϕ k = 2kπ/N is then equal to the maximum signal value in that direction, where the phase of the signal's derivative is ψ(t) = 2kπ/N + π/2. Thus, the mean for the first algorithm reads

m B1 (t) = 1 N k x ψ -1 2kπ N + π 2 . (6)
In the limit where the number of directions tends to infinity, this results in

m B1 (t) = 1 2π 2π 0 x(ψ -1 (ϕ))dϕ = 1 2π T 0 x(t) dψ dt dt, (7)
which is simply the mean of the signal over a period weighted by dψ(t)/dt > 0, where the weighting conveys the fact that the distribution of sampling points on the tube section is denser where the curvature is larger. Likewise, the same reasoning for the second algorithm results in

m B2 (t) = 1 π T 0 e iψ(t) Re e -iψ(t) x(t) dψ dt dt, = m B1 (t) + 1 2π T 0 e 2iψ(t) x * (t) dψ dt dt, (8)
and hence m

B2 (t) = m B1 (t) since T 0 e 2iψ(t) x * (t) dψ dt dt = i 2 -e 2iψ(t) x * (t) T 0 + T 0 dx * dt e 2iψ(t) dt = i 2 T 0 r(t)e iψ(t) dt = i 2 T 0 dx dt dt = 0 (9)
Thus, in this very simple case, the mean is in fact the same for both algorithms and therefore such a simple signal is a fixed point of both sifting operators iff the integral ( 7) is close to zero.

More generally, the outputs of the two proposed algorithms are very similar when the data clearly contains rotating components, but they may differ significantly when they fail to extract rotating components. Notice, though, that rarely one method succeeds in retrieving rotating components when the other one fails.

IV. ILLUSTRATION

A typical application of the proposed algorithms is proposed in Fig. 3. The data is a position record from an acoustically tracked, neutrally buoyant subsurface oceanographic float, one of a number deployed in the eastern subtropical North Atlantic Ocean in order to track the motion of dense salty water flowing out from the Mediterranean Sea during the "Eastern Basin" experiment [START_REF] Richardson | Tracking three Meddies with SOFAR floats[END_REF]. The data is available online from the World Ocean Circulation Experiment Subsurface Float Data Assembly Center (WFDAC) at http://wfdac.whoi.edu. Looping trajectories are indicative of intense swirling currents around an isolated packet of Mediterranean Sea water. Such structures, called "coherent vortices", are frequently observed in the ocean [START_REF] Mcwilliams | Submesoscale coherent vortices in the ocean[END_REF] and are more generally a ubiquitous feature of rotating turbulent fluids [START_REF]The vortices of geostrophic turbulence[END_REF].

Applied to such signals which a priori contain meaningful rotating components, the output of the bivariate extensions typically provide the given decomposition, where the rotations that were apparent in the original signal have been isolated in separate components. As mentioned earlier, though, not all components are rotating but primarily the first two ones, which correspond to the presumed coherent vortex. On this example however, we do not precisely know what information can be extracted from the decomposition yet, but we expect the large scale non-rotating components to reveal useful information regarding the background fluctuations determining the vortex position while the rotating components can a priori be used to extract finer informations, such as amplitude, angular frequency,etc. . . Such advanced study of the rotating components has already been performed using wavelet ridges to extract the coherent vortex signal [START_REF] Lilly | Wavelet ridge diagnosis of time-varying elliptical signals with application to an oceanic eddy[END_REF]. A comparative study of the EMD bivariate extensions for similar purposes is currently under investigation.

Fig. 1 .

 1 Fig. 1. The principle of the bivariate extensions. (a) A composite rotating signal. (b) The signal enclosed in its 3D envelope. The black thick lines stand for the envelope curves that are used to derive the mean. (c) Rapidly rotating component. (d) More slowly rotating component corresponding to the mean of the tube in (b).

Fig. 2 .

 2 Fig. 2. of the mean of the envelope for each algorithm. The accuracies of the estimation of the sampling points and of their mean are represented schematically as the thick lines of variable length under each point. (the provided script allows to easily test other configurations)


Valentin Genuer

Olivier Gal

Jérémy Méteau

Pierre Marcoux

Emmanuelle Schultz

Éric Lacot

Max Maurin

Jean-Marc Dinten

Optical elastic scattering for early label-free identification of clinical pathogens

Keywords: Optical elastic scattering, bacteria identification, bacterial colony, early screening, optical elastic backwardscattering, SVM, Fresnel diffraction patterns, MRSA

We report here on the ability of elastic light scattering in discriminating Gram+, Gram-and yeasts at an early stage of growth (6h). Our technique is non-invasive, low cost and does require neither skilled operators nor reagents. Therefore it is compatible with automation. It is based on the analysis of the scattering pattern (scatterogram) generated by a bacterial microcolony growing on agar, when placed in the path of a laser beam. Measurements are directly performed on closed Petri dishes.

The characteristic features of a given scatterogram are first computed by projecting the pattern onto the Zernike orthogonal basis. Then the obtained data are compared to a database so that machine learning can yield identification result. A 10-fold cross-validation was performed on a database over 8 species (15 strains, 1906 scatterograms), at 6h of incubation. It yielded a 94% correct classification rate between Gram+, Gram-and yeasts. Results can be improved by using a more relevant function basis for projections, such as Fourier-Bessel functions. A fully integrated instrument has been installed at the Grenoble hospital's laboratory of bacteriology and a validation campaign has been started for the early screening of SA and MRSA (Staphylococcus aureus, methicillin-resistant S. aureus) carriers.

Up to now, all the published studies about elastic scattering were performed in a forward mode, which is restricted to transparent media. However, in clinical diagnostics, most of media are opaque, such as blood-supplemented agar. That is why we propose a novel scheme capable of collecting back-scattered light which provides comparable results.

INTRODUCTION

Improved, fast and accurate detection and diagnostic technologies are needed in order to address present and future needs for rapid identification of contaminant or pathogens such as bacteria and yeasts for both industrial microbiological control and clinical diagnostics [START_REF] Newell | Food-borne diseases -the challenges of 20 years ago still persist while new ones continue to emerge[END_REF].

Although traditional culturing methods employing classical microbiological techniques are considered as gold standards because they provide high accuracy rates, they suffer from lengthy process due to the need of multiple preparation steps; taking up to 3-7 days to obtain negative results and additional 2-3 days to confirm a positive result by using biochemical or molecular techniques [START_REF] Bhunia | One Day to One Hour: How Quickly Can Foodborne Pathogens Be Detected?[END_REF]. Thus, since a decade, there is a growing interest for alternative methods based on light scattering which potentially have considerable advantage over the standards because they allow fast and non-invasive direct on agar analysis of bacterial colonies. The overall process of those techniques based on light scattering starts by noting that differences in refractive index, biomass, characteristics sizes and internal arrangement of bacteria inside the colonies leads to different forward scattering patterns when shooting it with a laser beam. Then with the help from learning algorithms it is possible to build a model based on a database containing numerous scattering patterns corresponding to different microorganisms and thus to provide discrimination results between them [START_REF] Bhunia | Light scattering, fiber opticand cellbased sensors for sensitive detection of foodborne pathogens[END_REF][START_REF] Banada | Label-free detection of multiple bacterial pathogens using light-scattering sensor[END_REF][START_REF] Suchwalko | Identification of bacteria species by using morphological and textural properties of bacterial colonies diffraction patterns[END_REF].

Advanced Bioimaging System was the first to launch in 2009 his BARDOT: a bacterial rapid detection system using optical light scattering developed at Purdue University which is able to provide precise identification of major foodborne pathogens. Around 90-99% for a discrimination at genus or species level for Listeria, Salmonella, Escherichia or Staphylococcus [START_REF] Bae | Label-free identification of bacterial microcolonies via elastic scattering[END_REF][START_REF] Tang | Light Scattering Sensor for Direct Identification of Colonies ofEscherichia coli Serogroups O26, O45, O103, O111, O121, O145 and O157[END_REF]. However his analysis throughput is limited by the incubation time needed to let the bacterial colonies reach a certain size, typically around 2 millimeters. Recently, a study by P.R. Marcoux et al reported some results about the development of a new kind of system using elastic forward scattering which provide good discrimination rates at species and strain level on bacterial and yeasts microcolonies at an early stage of growth. For example an 87% recognition rate for a strain level discrimination of four strains of Escherichia coli after 6h of incubation at 37°C was reached. The novel setup they proposed combined two optical paths in order to acquire images of microcolonies both in direct space (morphology insight) and in reciprocal space (diffraction pattern). Furthermore they used a specific laser beam profile to match the variability in size of the different microorganisms after 6h of growth (between 20µm and 250µm) [START_REF] Marcoux | Optical forward-scattering for identification of bacteria within microcol onies[END_REF].

The present study inspired by these recent developments reports on the improvements we brought both to the optical design and the image analysis so as to increase classification performances and also to address clinical applications involving non-transparent growth culture media.

It first deals with the optical design and the results obtained in discriminating Gram+, Gram-and yeasts at 6 hours with our forward-scattering system; the screening of SA and MRSA carriers will also be explained. Then it details the improvement that can be achieved by working on the analysis process, especially by changing the orthogonal basis on which the images are projected to build the features vectors.

Yet, the technique is still restricted to transparent culture media due to transmission acquisition geometry, and is consequently not adapted to coloured media such as blood-supplemented agar media. Those opaque media are of common use in clinical microbiology and up to now, there is not any system based on elastic scattering that can provide comprehensive scattering patterns from microcolonies growing on such culture media. A study, carried out in 2004 by S. Guo showed that forward scattering was more promising than back-scattering method because they could not build a system capable of providing patterns of good quality [START_REF] Guo | Optical scattering for bacterial colony detection and characterization[END_REF]. This is why we propose a novel system that could enlarge the application domain of elastic light scattering methods to the use of opaque culture media by providing an optical architecture capable of overcoming this challenge and thus collecting comprehensive patterns from light back-scattered by microcolonies directly on their commercial Petri dishes. We show here the solutions we provide to build such a system and we also show the results obtained after testing its discrimination performance at the species or strain level on different bacterial and yeasts species. For the testing, a database was acquired using a common opaque culture media (blood-supplemented agar COS) used in clinical microbiology to make fastidious organisms grow.

MATERIALS AND METHODS

Sample preparation

Bacterial strains were obtained from KwikStik (Microbiologics, St. Cloud, MN) lyophilized ATCC reference strains. Starting from a culture on COS or TSA (bioMérieux) after 24h of incubation (37°C), a 0.5 McF suspension (Densimat, bioMérieux) was prepared in API NaCl 0.85 % Medium (bioMérieux). This 0.5 McF was then immediately diluted at 1/120 in NaCl 0.85 %. 50µL of this diluted suspension were finally plated on agar medium (90mm petri dish, bioMérieux), either on ChromID S. aureus or TSA, for 6h of additional incubation and scattering experiments.

Instruments

The scatterograms were acquired using the forward-scattering setup depicted in Figure 1. A 532 nm laser (Spectra Physics Excelsior 532-50-CDRH) produces a coherent collimated light beam with a power decreased to less than 1mW thanks to an optical density. A pair of achromatic lenses is used to give this beam a particular shape, allowing us to reach a range of laser spot sizes from 10µm to 300µm by moving the sample along the beam path. A first video camera (Pixelfly QE) is used to acquire the forward scatterograms resulting from the interaction between the laser beam and the targeted microcolony. A second video camera coupled with achromatic lenses is used to get a microscope-like image of the sample. This bright field image is obtained when the sample is put at the focus point of the 2-lenses-system and illuminated by a ring LED placed above it. This configuration is convenient because it provides both a direct image of the targeted microcolony and its corresponding forward scatterogram. The alignment of the sample with the laser beam, and the adjustment of the laser spot size on the sample are achieved using horizontal and vertical translation stages (PI Micos VT75).

Recently, this setup has been integrated in a transportable instrument, the so-called MICRODIFF system, depicted in Figure 2. It is currently being tested at the Grenoble hospital's laboratory of bacteriology for the screening of the SA and MRSA (S. aureus and methicillin-resistant S. aureus) carriers. Since the Staphylococci that were used for the database grow slowly, we had to replace the lens above the sample by a microscope objective (10x, NA: 0.25) so as to get a laser spot size that could reach the microcolonies sizes after a 6h incubation time. At this time of growth, indeed, some strains were only around 4 µm long.

To enlarge the application domain of elastic scattering to clinical diagnostic purpose needing Petri dishes filled with opaque culture medium, we build another system, depicted in Figure 3, able to acquire scatterograms backward reflected by the sample. The same laser source coupled with achromatic lenses is used to produce a light beam with a particular shape; the same principle is used to vary the laser spot size on the sample: we move the sample across the beam path. There is a direct imaging system too composed of a CMOS video camera and a couple of achromatic lenses (orange path in Figure 3). The focus is obtained by placing the sample at the focal plane of the lens near it. To acquire scatterograms backward reflected we used a quarter-wave plate coupled with a polarizing beam splitter (blue path in Figure 3). The light issued from the laser source is linearly polarized and is first reflected toward the sample thanks to the polarizing beam splitter (Thorlabs, CM1-PBS251). As it passes through the quarter-wave plate, it becomes rotationally polarized. The reflection on the sample inverts the direction of polarization and as the light passes back through the quarter-wave plate it becomes linearly polarized but perpendicularly to the initial one. Thus it is now transmitted through the polarizing beam splitter and then collected by an achromatic lens and finally acquired on a CMOS video camera. This system could reach a spot size range from 30 µm to 100 µm and had a numerical aperture around 0.15. This setup has also recently been upgraded by adding a microscope objective (20x magnification, NA: 0.45) to increase the numerical aperture and thus get more diffraction rings and also get a microscope image with a greater magnification.

The main drawback of adding optical elements between the polarizing beam splitter and the sample is getting specular reflections that saturate the CMOS sensor. This is managed by tilting some elements, such as the quarter-wave plate while adjusting the polarization with the quarter-wave plate.

Image Analysis

Zernike moments invariants

As a first choice for features-extraction from images, the use of Zernike moments invariants was quite relevant because of the typically circularly shaped scattering patterns. Basically, those scatterograms were initially centered in the middle and histograms were equalized. Then Zernike moments invariants were calculated up to the 20 th order.

Zernike moments are computed from a complete orthogonal basis set of polynomials defined on the unit disk [START_REF] Teague | Image analysis via the general theory of moments[END_REF]. The expression of these Zernike orthogonal basis functions in polar coordinates is classically:

𝑉 𝑛𝑚 (𝑟, 𝜃) = 𝑅 𝑛𝑚 (𝑟)exp (𝑗𝑛𝜃),
where (r,θ) are defined on the unit disk, n is a non-negative integer, |m| is less or equal to n and n-|m| must be even. R nm is the Zernike radial polynomial, defined in detail in [START_REF] Khotanzad | Invariant image recognition by Zernike moments[END_REF] for example.

Zernike moments are calculated after projecting the images on this orthogonal basis. Thus if we call f the digital image function, the Zernike moment of radial order n with azimuthal repetition m is expressed as:

𝑍 𝑛𝑚 (𝑟, 𝜃) = 𝑛+1 𝜋 ∑ ∑ 𝑓(𝑥, 𝑦)𝑉 𝑛𝑚 * (𝑟, 𝜃)𝑑𝑥𝑑𝑦 𝑦 𝑥 , 𝑥 2 + 𝑦 2 ≤ 1.
As rotation-invariant features we took the magnitudes of Zernike moments as they remain unchanged under rotation. Computing these moments up to the 20 th order gave us features-vector of 121 components per image which then were used to train a supervised machine learning algorithm. The main idea was to build a model that correctly mapped the features-vector to the corresponding sample (e.g. bacterial strain). The process, previously detailed [START_REF] Marcoux | Optical forward-scattering for identification of bacteria within microcol onies[END_REF], follows three main steps: first the acquisition of a large database composed of hundreds of images per class; then the training of the learning algorithm, using an SMO (Sequential Minimal Optimization) which is a particular algorithm used to train SVM (Support Vector Machines) while reducing calculation time; finally a cross-validation to test the eventual classifier errors and summarize them in a confusion matrix. The WEKA data mining software was used for this last task.

A first observation about Zernike radial polynomials properties is that at the order n we can at most obtain n/4 complete oscillations. This limits at n/4 the number of rings in scatterograms that can be caught by the projection. For example (Figure 4a), at the order n=20 we can hope to retrieve at most 5 rings. This is too few compared to the actual number of rings present in the scatterograms we observed, and increasing the order would have drastically increased the features-vector size. Moreover we can note that intensity always reaches his maximum for r=1 (i.e. at the unit circle boundaries) while we observed on our images that intensity was rapidly decreasing with r. Thus it could be problematic in terms of expansion convergence. This is why we decided to experiment another projection basis.

Fourier-Bessel moments invariants

Fourier-Bessel moments are computed from a basis of functions based on the Bessel function of the first kind. The functions:

𝐽 𝑛 (𝛼 𝑛𝑚 𝑟)𝑒 ±𝑗𝑛𝜃 ,
where n ≥ 0, m ≥ 1 and α nm is the m th root of the J n Bessel function, form an orthogonal basis on the unit disk r ≤ 1 [START_REF] Zhao | Fourier-Bessel rotational invariant eigenimages[END_REF]. Calling f the image function, Fourier-Bessel moments are defined in polar coordinates as follows:

𝐹 𝑛𝑚 = 1 2𝜋𝑐 𝑚 ∫ ∫ 𝑓(𝑟, 𝜃)𝐽 𝑛 (𝛼 𝑛𝑚 𝑟)𝑒 -𝑗𝑛𝜃 𝑟𝑑𝑟𝑑𝜃 1 0 2𝜋 0 ,
where m ≥ 1, n = 0, ±1, ±2, … is the moment order, c m is a normalization constant and α nm is the m th root of the J n Bessel function. As we did with the Zernike moments, we took the magnitudes of the Fourier-Bessel moments to build our rotation invariants.

It was previously shown [START_REF] Xiao | Image analysis by Bessel-Fourier moments[END_REF] that Fourier-Bessel moments take advantages of the zeros repartition of the Fourier-Bessel radial function J n (α nm r), which are uniformly distributed over the interval 0 ≤ r ≤ 1 (Figure 4b for n = 20) unlike the zeros of the Zernike radial polynomials which are mainly located at the unit circle boundaries (Figure 4a).

RESULTS AND DISCUSSION

Gram+/Gram-/Yeasts discrimination

Figure 5 presents some scatterograms of the first database acquired in this work. 1906 scatterograms were acquired over the course of 3 months on Gram+, Gram-and yeasts (more than 120 scatterograms per strain) for a total of 15 different strains cultured in the conditions described above. Petri dishes lids remained closed during the whole acquisition process. The features-vectors were computed using 121 Zernike moments and the SMO was trained. The 10fold cross-validation showed a global correct classification rate of 75% down to the strain level (details not shown). A more relevant classification rate of around 95% was obtained for the discrimination between Gram+, Gram-and yeasts. The confusion matrix is detailed in Table 1.

Early screening of SA/MRSA carriers at 6h

A validation campaign has been started at the Grenoble's hospital laboratory of bacteriology on the early screening of Staphylococcus aureus (SA) and methicillin-resistant S. aureus (MRSA) carriers. The acquisitions were made at 6h of growth with the fully integrated instrument version of the forward-scattering system MICRODIFF. A specific plate, the ChromID S. aureus was used. In fact this is a bi-plate filled with a specific medium for the growth of Staphylococci. One side contains a growth medium enriched with methicillin to detect resistant Staphylococci and the other side is only filled with the growth medium. The database, which is still in acquisition, is restrained to clinical relevant Staphylococci growing on ChromID and some interfering strains (i.e. non-Staphylococci strains which are growing on ChromID, e.g. some gram-like Stenotrophomonas matlophilia and yeast like Candida tropicalis). At this time we are able to get an average correct classification rate of 75.5% ± 0.2% (average on ten 10-fold cross-validations) between S. aureus and other Staphylococci (using 121 Zernike moments). This is a promising result at 6h of growth. Nevertheless this clinical application shows the difficulty of the exercise in terms of bacterial microcolony size. As we can see in Figure 6, scatterograms acquired on Staphylococci after only 6h of growth are quite similar and simple. This is principally due to the small size of the microcolonies and also the low biomass, since those strains are growing slowly. Average correct classification rate: 87.0% ± 0.6%* * standard deviation over ten 10-fold cross-validations Table 3 Confusion matrix on E. coli strain-level discrimination at 6h of growth (classification rates in %), using 500 Fourier-Bessel moments. Same database as for Table 2. Average correct classification rate: 93.0% ± 0.7%* * standard deviation over ten 10-fold cross-validations

Figure 1

 1 Figure 1 Schematic of the forward-scattering experimental setup for the acquisition of scatterograms on transparent growth media. Abbreviations: L: achromatic lens, BS: beam splitter, RL: ring LED. CMOS sensors: (a) for the direct space imaging, (b) for the acquisition of the scatterograms.

Figure 2

 2 Figure 2 Photograph of the integrated MICRODIFF instrument for fast identification of clinical pathogens using elastic light scattering. Dimensions: 63x21x21 cm3, 10kg.

Figure 3

 3 Figure 3 Schematic of the backward-scattering system. Abbreviations: L: achromatic lens, BS: beam splitter, PBS: polarizing beam splitter, RL: ring LED. CMOS sensors: (b) for the direct space imaging, (a) for the acquisition of the scatterograms.

Figure 4 (Figure 5

 45 Figure 4 (a) Zernike radial polynomials R nm (r) for m=10. (b) Radial functions J 1 (α n r) for the Fourier-Bessel projections with azimuthal order 1. We can note the more convenient zeros distribution of the Fourier-Bessel radial functions.

Figure 6

 6 Figure 6 Scatterograms of microcolonies acquired after 6h of growth at 37°C (on ChromID S. aureus) with the forward scattering prototype at the Grenoble's hospital laboratory of bacteriology during the validation campaign on the early screening of SA/MRSA carriers.

Figure 7

 7 Figure 7 Images of an E. coli GFP microcolony acquired on the backward-scattering setup at 6h of growth on LB + ampicillin. (a) Direct space imaging. (b) Direct space imaging plus laser spot targeting the micro-colony. (c) Scatterogram acquired on the setup. (d) Microscope image with a 50x objective, AxioVision Zeiss.

Figure 8

 8 Figure 8 Comparison of the image reconstructions obtained using Zernike and Fourier-Bessel expansions. (a) Scatterogram of an E. coli ATCC 25922 microcolony acquired with the forward-scattering system MICRODIFF after 6h at 37°C. (b) Reconstruction after projection on the 121 first Zernike moments (maximal radial order: 20, maximal azimuthal order: 20).Quadratic error: 66.3%. (c) Reconstruction from the projection on 120 Fourier-Bessel moments (maximal radial order: 24, maximal azimuthal order: 4). Quadratic error: 37.7%

Figure 9

 9 Figure 9 Reconstruction error decrease with increasing number of moments, for Zernike (blue *) and Fourier-Bessel (red o) expansions.

Table 1

 1 Confusion matrix on Gram+/Gram-/Yeasts discrimination at 6h of growth (classification rates in %), using 121 Zernike moments. Database of approximately 1900 scattograms from 15 strains. Average correct classification rate: 94.7% ± 0.1%* *standard deviation over ten 10-fold cross-validations

			Classified as
		Gram+	Gram-	Yeasts
	Gram+	92.7	3.4	3.9
	Gram-	2.6	92.4	5.0
	Yeasts	0.4	2.8	96.8

Table 2

 2 Confusion matrix on E. coli strain-level discrimination at 6h of growth (classification rates in %), using 121 Zernike moments. Database of 400 scatterograms from 4 strains.

			Classified as
		ATCC 25922 ATCC 35421 ATCC 11775 ATCC 8739
	ATCC 25922	96	1	1	2
	ATCC 35421	7	80	4	9
	ATCC 11775	6	6	80	8
	ATCC 8739	5	4	2	89

Validation of the backward scattering system on four E. coli strains

We decided to validate our backward scattering system by acquiring a database on four Escherichia coli strains commonly associated with human infections (ATCC 25922, 35421, 11775 and 8739). We used the same four strains as P. Marcoux et al [START_REF] Marcoux | Optical forward-scattering for identification of bacteria within microcol onies[END_REF] on their forward scattering system to test the strain-level discrimination. Commercial COS Petri dishes (Columbia agar +5% sheep blood) were inoculated and incubated during 6h at 37°C. Then the acquisition was performed using the instrument described above. The direct space imaging unit for fast colony localization and targeting combined with short acquisition time (around 10 ms) allowed to acquire around 100 scatterograms per strain in a week, for a database of 400 images. Figure 7 shows the scatterogram acquired in reflection on a bacterial microcolony. It also displays the microscope image of the microcolony compared to the one observed with the bright field unit of the setup.

To rigorously compare the performances of the backward scattering system to the forward scattering system ones, we used the same features-vectors (121 Zernike invariants moments) to describe the images. Then the same SMO learning algorithm coupled with a 10-fold cross-validation was used. The confusion matrix showed an overall correct classification rate of 87% (Table 2), which is close to the 90% obtained with the forward scattering system and thus validates the new instrument. Promising results were also obtained on the discrimination of three Candida yeasts species (ATCC 14053, 2001 and 18303) where we could reach an overall classification rate of 94% (details not shown).

Analysis improvements with Fourier-Bessel moments

We made two main observations about the use of Fourier-Bessel moments: first we looked at the accuracy of the projection by comparing Zernike and Fourier-Bessel projections on different scatterograms and then we tried to evaluate how the number of moments used for the features-extraction influenced the classification rates. Figure 8 shows the comparison between the projections of a scatterogram on Zernike basis and on Fourier-Bessel basis. Once the projections were made, we computed the reconstruction and calculate the relative error (to the original image). Approximately the same amount of coefficients was used: 121 for the Zernike projection (maximal radial order: 20, maximal azimuthal order: 20) and 120 for the Fourier-Bessel one (maximal radial order: 24, maximal azimuthal order: 4). It showed a better reconstruction accuracy for Fourier-Bessel, 37%, against 66% for Zernike with approximately the same number of moments.

We also showed that increasing the number of moments gives a more faithful reconstruction and thus a lower reconstruction error. Figure 10 illustrates that at equal number of moments, using Zernike basis leads to higher reconstruction error rates. On this example we reached a plateau at around 30% reconstruction error with the Fourier-Bessel basis. This level was reached from 160 Fourier-Bessel moments. It could be due to the noise level and the presence of local defects in the image. The Zernike reconstruction error seemed to decrease steadily but was always far higher than the Fourier-Bessel.

We used the Fourier-Bessel moments on the 4-strains E-coli database acquired with the backward system. We reached an average classification rate of 93.0% ± 0.7% (average over ten 10-fold cross-validations; confusion matrix given in Table 3). Compared to the 87% obtained when using Zernike moments, it corresponds to a decrease of the classification error rate by a factor 1.9.

CONCLUSION

The MICRODIFF instrument was employed to perform fast identification of bacterial pathogens down to the strain level. Several scatterograms databases were acquired at an early stage of growth in a forward or backward scattering mode. A correct classification rate of 94% was achieved on the Gram+/Gram-/Yeasts discrimination using the forward scattering system and a SVM classifier in a 10-fold cross-validation.

The backward scattering system allowed us to reach 87% on the classification of 4 E. coli strains. This has validated the optical reflection design for the identification of pathogens directly on opaque agar media. A novel setup is now being developed to combine both forward and backward scattering acquisitions and also the direct bright field imaging unit.

We also showed we could improve these results by using Fourier-Bessel moments instead of Zernike moments for the features extraction from scattering patterns. The correct classification on the 4 E. coli strains database rate was increased to 93%, corresponding to a reduction of the classification error rate by a factor 1.9. This shows that decreasing the reconstruction error and thus finding the best way to faithfully describe the scatterograms leads to higher classification rates. Further work is needed to more fully understand how much we can still gain in classification performance and also to determine whether the classification process relies on the whole features vector or only on few components.


Dongzhe Li

Cyrille Barreteau

Martin R Castell

Fabien Silly

Alexander Smogunov

Out-versus in-plane magnetic anisotropy of free Fe and Co nanocrystals: tight-binding and first-principles studies

Keywords: numbers: 71.15.-m, 75.10.Lp, 75.50.Ss, 75.70.-i, 75.75.Lf, 68.47.Jn

) is out-of-plane, it is in-plane for Co(001). This has direct consequences on the magnetic reversal mechanism of the nanocrystals. Indeed, the very high uniaxial anisotropy of Fe nanocrystals makes them a much better potential candidate for magnetic storage devices.

I. Introduction

Last decades, higher storage densities were achieved by reducing the magnetic grains down to nanoscale. However, the magnetic stability of a nano-object decreases proportionally to its size and the ultimate limit is reached when the thermal fluctuation overcomes the energy barrier to swicth the global magnetization of the system. The most crucial issue in exploring ultimate density data storage (e.g., highdensity magnetic recording [START_REF] Albrecht | [END_REF] or spintronic devices) is magnetic anisotropy energy, which is defined as the change of total energy associated to a change in the direction of magnetization. One of the challenge in this route towards high magnetic density storage is evidently to be able to synthetize well ordered arrays of magnetic nanocrystals with as large magnetization and magneto-crystalline anisotropy as possible. The magnetic anisotropy energy of magnetic nanocrystals (e.g., Fe, Co and Ni etc) is indeed a subject of intense study both experimentally [2][3][4][5] and theoretically [6][7][8][9][10][11] but the ability to grow well defined magnetic crytalline nanostructures is also a major issue [12][13][14][15][16][17][18] . This is especially the case for Fe and Co nanostructures, that can adopt various crystalline bulk structures, in particular the body-centered cubic (bcc) and face-centered cubic structure (fcc) structure in low dimensions 15,16 . The closepacked and lowest-energy facet for bcc structure is the (110) facet whereas it is the (111) facet for the fcc structure. This is the reason for the (110) facets appearing in bcc nanocrystals (in fact for Fe the surface energies of (001) and (110) orientations are almost the same) and the (111) facets appearing in fcc nanocrystals. The nanocrystal magnetic properties will therefore not only depend on the bulk structure but also on the facet structure and area.

The magnetic anisotropy contains two different parts: the first part is long-range magnetic dipole-dipole interaction which leads to so-called shape anisotropy, while the second one is referred to magnetocrystalline anisotropy energy (MAE) originates from the spin-orbit coupling (SOC) 19 . The latter effect is a quantum effect (of relativistic nature) that breaks the rotational invariance with respect to the spin quantization axis. Therefore, if SOC is included, the energy of the system depends on the orientation of the spin with respect to the crystallographic axis.

The value of MAE per atom, is extremely small in bulk (some µeV), but can get much larger in nanostructures 20,21 (some meV) due to reduced dimensionality. From the point of view of theory, there are two different methods which are used extensively in the literatures for MAE calculations: i) fully relativistic self-consistent field (SCF) calculations, ii) Force Theorem (FT) [22][23][24] . Assessing the MAE for systems containing hundreds of atoms by the former approach is especially challenging, since it requires a well-converged charge density as well as a consuming computational time. In the latter method, the MAE is given by the band energy difference (instead of total energy difference) obtained after a one-step diagonalization of the full Hamiltonian including SOC, starting from the self-consistent scalar relativistic density/potential. This approach is not only computationally efficient but also numerically very stable since the self-consistent effect with SOC could be ignored.

The objective of this paper is to investigate the MAE of Fe and Co nanocrystals, that can be grown experimentally by epitaxy, using tight-binding (TB) as well as first principles calculations in the Density Functional Theory (DFT) framework. The nanocrystals adopt a truncated-pyramid shape on a reconstructed SrTiO 3 (001) substrate but have however a different bulk structure (bcc and fcc). In our theoretical study a particular emphasis is devoted to the local analysis of MAE. In particular, it has been found that the main contribution to the MAE comes from the basal (001) facet of pyramids. This results in strong out-of-plane (in-plane) anisotropy for Fe (Co) nanocrystals, in agreement with the study of thick Fe(001) and Co(001) slabs.

The paper is organized as follows. In Sec. II we present the experimental and theoretical methods used in this work. In Sec. III, we first present Scanning Tunneling Microscope (STM) observation of Co nanocrystals on SrTiO 3 (001) substrate and illustrate the results of TB and DFT calculations for arXiv:1407.5830v1 [cond-mat.mtrl-sci] 22 Jul 2014 Co(001) and Co(111) slabs. After that, the MAE of free Fe and Co nanocrystals will be discussed. Finally, the conclusions will be presented in Sec. IV.

II. Methodolody

In the following sections we will first briefly present the main ingredients of the experimental set-up to grow the cobalt nanocrystals on the SrTiO 3 (001) substrate. Then in the second part the theoretical model to calculate the MAE will be presented starting by the TB Hamiltonian and then the DFT approach. In the case of the DFT formalism we will essentially concentrate on the implementation of the Force Theorem which we have incorporated in the Quantum-ESPRESSO (QE) 25 package.

A. Experimental we use SrTiO 3 (001) crystals doped with 0.5% (weight) Nb. The crystals were epi-polished (001) and supplied by PI-KEM, Surrey, UK. We deposited Co from an e-beam evaporator (Oxford Applied Research EGN4) using 99.95% pure Co rods supplied by Goodfellow, UK. Our STM is manufactured by JEOL (JSTM 4500s) and operates in ultra high vacuum (10 -8 Pa). We used etched W tips for imaging the samples at room temperature with a bias voltage applied to the sample. SrTiO 3 (001)-c(4×2) was obtained after Ar + bombardment and annealing in UHV at 600 • C for 2 hours. STM images were processed and analyzed using the home made FabViewer application 26 .

B. Theoretical

Magnetic tight-binding model

In this section, we briefly describe our magnetic tightbinding model (more details can be found in our previous publications 27,28 ). The hamiltonian is written as follows:

H = H TB + H LCN + H Stoner + H SOC (1)
Where H TB is a standard "non-magnetic" TB hamiltonian which form is very similar to the one introduced by Mehl and Papaconstantopoulos 29 , H LCN is the term ensuring a local charge neutrality, H Stoner is the Stoner-like contribution that controls the spin magnetization and H SOC corresponds to spin-orbit coupling that operates on d orbitals only.

The total energy should be corrected by a double counting term due to inter-electronic interactions introduced by local charge neutrality and Stoner interaction, explicitly:

E tot = E b -E dc = α f α α - U 2 i [N 2 i -(N 0 i ) 2 ]+ 1 4 i,λ I λ M 2 iλ , (2)
where E b = α f α α is the band energy, f α is the Fermi-Dirac occupation of state α, N i and M i are the charge and the spin moment of site i, respectively. N 0 i is the valence charge, U is the parameter imposing the local charge neutrality, and I λ is the Stoner parameter of the orbital λ (λ = s, p, d).

All the parameters of TB hamiltonian are fitted on bulk ab initio data: bandstructure, total energy, magnetic moment etc. The value of the Stoner parameter I d is taken equal to 0.88 meV for Fe and 1.10 meV for Co. The spin-orbit constant ξ d is also determined by comparison with ab initio bandstructure and we found that 60 meV and 80 meV are very good estimates for Fe and Co, respectively.

The MAE, in a very good approximation, is calculated by using the Force Theorem (FT) [22][23][24] : first, a self-consistent field (SCF) collinear calculation without SOC is done followed by the rotation of the density matrix in the right spin direction; next, a non-SCF non-collinear calculation with SOC is performed. The MAE is obtained as the difference of band energies, E 1

b -E 2 b , between two spin moment directions, 1 and 2. The correct decomposition of total MAE over different atomic sites i can be done within so-called "grand canonical" formulation 24 :

M AE i = E F (E -E F )∆n i (E)dE (3)
where

∆n i (E) = n 1 i (E) -n 2 i (E)
is the change in the density of states at atom i for different spin moment orientations. The Fermi energy E F of SCF calculation without SOC must be subtracted from all energies in order to suppress the trivial contribution to the local MAE due to charge redistribution as discussed in Ref. 24 .

Density Functional Theory (DFT) calculations

We perform ab initio DFT calculations using the planewave electronic structure package Quantum-ESPRESSO (QE). 25 The spin-orbit coupling (SOC), crucial for magnetocrystalline anisotropy, is taken into account via fullyrelativistic pseudo-potentials (FR-PPs) 30 , describing the interaction of valence electrons with ions, which are in turn generated by solving atomic Dirac equations for each atomic type. We have implemented the Force Theorem in QE in the same two-step way as described above for TB model: i) SCF calculation with scalar-relativistic PPs (without SOC) is performed to obtain the charge density and the spin moment distributions in real space; ii) spin moment is globally rotated to a certain direction followed by a non-SCF calculation with FR-PPs (with SOC). The change of band energy between two spin moment directions gives, as above, the total MAE. The total MAE is decomposed over different atoms i in the slightly different way:

M AE i = E 1 F (E -E 2 F )n 1 i (E)dE - E 2 F (E -E 2 F )n 2 i (E)dE, (4)
where the Fermi level of one of magnetic configurations (we have chosen the second one), E 2 F , is substracted under integrals and exact Fermi levels for two configurations are used as the limits of integration. This way we avoid the reference to electronic levels of a system without SOC, since the PPs with and without SOC are not generally correlated and can produce an arbitrary shift of levels. Due to total charge conservation in this "canonical" approach, the sum of M AE i over all atoms gives exactly the total MAE while for the "grand canonical" scheme, Eq. 3, it was, in principle, only approximate. The descrepancy between "grand canonical" and "canonical" formulations within TB approach is, however, very tiny since the effect of SOC on the Fermi level is negligable in the case of Fe or Co composed materials. Since QE gives an access to real space wave-functions it is natural to define also the space-resolved MAE as:

M AE(r) = E 1 F (E -E 2 F )n 1 (r, E)dE- E 2 F (E -E 2 F )n 2 (r, E)dE, (5)
where the LDOS is computed via electron wave-functions in the usual way,

n 1,2 (r, E) = α |Ψ 1,2 α (r)| 2 δ(E -ε 1,2 α ).
Once again, the integral of MAE(r) over all the space will give exactly the total MAE.

III. Results and discussion

In this section we will first briefly present the structural characterization of supported Co nanocrystals using STM. We then present the results of our calculations on slabs of fcc Co with orientations (001) and (111) corresponding to the facets of the nanocrystals. Next, in the second part we consider Fe and Co nanocsrystals in form of truncated pyramids with the same length to height ratio as in the experiments.

A. STM Observations

The SrTiO 3 (001)-c(4×2) surface 31 is used for cobalt deposition. The c(4×2) reconstruction was verified by STM and LEED before deposition. Fig. 1 shows the topography of the SrTiO 3 (001)-c(4×2) surface following deposition of 3 monolayers (ML) of Co on a substrate heated to 320 • C followed by a subsequent 50 minute anneal at 350 • C. The Co has selfassembled into similarly sized nanocrystals. Cobalt usually adopt a hcp bulk structure but STM shows that Co nanocrystals have a fcc structure with the shape of a truncated pyramid (a square top surface and a square base). The Co nanocrystals have. The side facets of the nanocrystals were measured at an angle of ∼54 • with respect to the substrate. This shows that Co is cubic packed and the nanocrystals have a (001) top facet and four (111) side facets. The interface is therefore a (001) plane and the interface crystallography is (001) Co (001) SrTiO3 , [100] Co [100] SrTiO3 . As a guide to the eye we have shown in Fig. 1 (inset) a schematic illustration of a truncated pyramid. The ratio of the length ( ) of the top square to the height (h) of the truncated pyramids as a function of volume is shown in Fig. 1b. The constant ratio of /h=1.48±0.13 suggests that these pyramidal nanocrystals have reached their equilibrium shape. The error in the ratio denotes the standard deviation of the measurements.

B. Calculations

As has been discussed above, fcc Co nanocrystals (Fig. 1) as well as bcc Fe nanocrystals 16 can be epitaxially grown on SrTiO 3 (001) substrate with a remarkable control of size, shape and structure. These crystals can contain up to several hundreds of atoms and have the form of truncated pyramids, as shown in Fig. 2, with a rather constant length-to-height ratio, l/h. They however adopt different bulk structure, i.e. the nanocrystal facets will therefore be different because the close-packed and lowest-energy facet for bcc structure is the (110) facet whereas it is the (111) facet for the fcc structure. It is expected that the MAE of such pyramids will be dominated by the surface composed of (001) and ( 110) or ( 001) and (111) facets for Fe and Co nanocrystals, respectively. It is therefore essential to estimate first the MAE of the bulk slabs of these orientations. We present below the results for fcc Co (001) and (111) slabs while similar results for bcc Fe slabs have already been reported recently (Ref. 24 ). 001) and (111) for Co, respectively. Their possible size and shape is controlled by length-to-height ratio, l/h, kept to ∼ 1.0 (Fe) and 1.41 (Co) which are close to experimental values, ∼ 1.20 (Fe) 16 and ∼ 1.48 (Co). The z axis was chosen to be normal to the pyramid base and the spin moment is rotated in the xz plane forming the angle θ with the z axis.

MAE of Co fcc (001) and (111) slabs

The Co slabs were constructed from fcc Co with a lattice parameter of a 0 = 3.531 Å found from ab initio calculations (which is close to the experimental value of a 0 = 3.548 Å ) and no atomic relaxations were performed. Fig. 3 shows thickness dependence of the total MAE of N-atom fcc Co slabs of (001) and (111) orientations. The results of both TB (N = 1∼20) as well as ab initio (N = 1∼10) calculations are presented. Note that the total MAE is obtained as total energy difference for M perpendicular or parallel to the atomic slabs. Explicitely, ∆E = E ⊥ tot -E tot . In the TB model, a mesh of 50 × 50 in-plane k points has been used for SCF calculations without SOC whereas the mesh was incresed to 70 × 70 in non-SCF calculations with SOC in order to provide a precision below 10 -5 eV. A Marzari-Vanderbilt broadening scheme with smearing parameter of 50 meV has been used. Ab initio DFT calculations were carried out with Quantum ESPRESSO package 25 using generalized gradient approximation (GGA) for exchangecorrelation potential in the Perdew, Burke, and Ernzerhof parametrization 32 . Full self-consistent calculations were performed with relativistic ultrasoft pseudo-potentials (no FT was employed here) and cut-off energies were set to 30 and 300 Ry for wave-functions and charge density, respectively. The mesh of 40 × 40 k points was used and the same smearing parameter and technique were employed.

We find a relatively good overall agreement between TB and DFT calculations. MAE oscillations for both slabs can be clearly seen even for quite thick slabs (similar results were recently reported for bcc Fe slabs 24 ). This kind of longrange oscillating behavior has been recently reported by experiments in thin ferromagnetic films (Fe and Co), and was interpreted in terms of spin-polarized quantum well sates. 33,34 . We notice further that for Co(001) slabs both calculations give rather similar results: the total MAE clearly favours in-plane magnetization with anisotropy energy around 0.6 meV/cell. In the case of Co(111), the MAE oscillates around zero in TB model while the DFT calculations predict rather small (compared to the (001) case) out-of-plane magnetic anisotropy. Note that our results compare rather well with DFT calculations in Ref. 35 done with LDA approximation for exchangecorrelation functional. We further study the local decomposition of MAE of ( 001) and ( 111) Co slabs made of 20 atomic layers (Fig. 4). Here, we used the FT in TB as well as in DFT approaches as described in the previous section. A qualitatively good agreement between TB and DFT calculations is again found for both slabs with the main discrepancy appearing for the surface layers, which indicates that the TB model is presumably less accurate for low coordinated atoms. Interestingly, for both (001) and ( 111) slabs these surface layers possess in-plane anisotropy. The local MAE site decomposition then shows damped oscillations converging towards a tiny bulk value. However, while the MAE of the (001) slab is strongly dominated by the outermost surfaces layer, this is not the case for the (111) slab where sub-surface layers cancel (and even overcome in the DFT case) the surface contribution. This leads to the large in-plane and rather small out-of-plane overall MAE for the (001) and (111) slabs, respectively, as it is reported in Fig. 3.

Free Fe and Co nanocrystals

The length-to-height ratio of different size of Fe and Co nanocrystals can be written

l/h = [2(n -1)]/(N -n) and l/h = (n -1)/[ √ 2(N -n)],
where N × N and n × n are the number of atoms in the first (bottom) and last (up) layers of the truncated pyramids. We then selected different sizes of Fe and Co nanocrystals with the length-to-height ratio of ∼ 1.0 (l/h = 1.0 for N =29, 135; 1.20 for N = 271; 1.14 for N = 620) and 1.41 close to the experimental value of 1.20±0.12 16 and 1.48±0.13, respectively. Since the MAE in the xy plane was found to be extremely small, we kept the magnetization always in the xz plane making the angle θ with the z axis. The MAE is defined as the change in the band energy between magnetic solutions with magnetization along the z and x axis, M AE = E z -E x . In Fig. 5, we plot the total MAE of Fe and Co nanocrystals of growing size calculated with TB approach. Different sign of MAE means that out-of-plane magnetization is favored in Fe while in Co the spin moment will be rotating in the easy xy plane which makes thus Fe nanocrystals a better candidate for magnetic storage applications. These results can One can see that the total MAE mainly originates from the lower (001) facet and its perimeter composed of least coordinated atoms. Therefore, in agreement with the previous analysis of (001) Co and Fe 24 slabs, this would favor the outof-plane/in-plane anisotropy for Fe/Co nanocrystals, respectively. We notice, moreover, that since nanocrystals of Co are much flatter then those of Fe (as Fig. 2 illustrates), which is a consequence of bigger length-to-height ratio for Co, in the case of Co nanocrystals also the upper (001) facet, containing more atoms, gives noticeable contribution to the overall MAE. We have also checked the total MAE in the xy plane but have found it extremely small, of amplitude about 3 meV and 0.8 meV for Fe (N=620) and Co (N=615) nanocrystals, respectively. As mentioned in Sec. I, another important contribution to magnetic anisotropy is the so-called shape anistropy energy. We have calculated it for biggest Fe (N=620) and Co (N=615) nanocrystals and have found rather small values, of about 5 meV and 2 meV for Fe and Co, respectively. Note that for both pyramids, the shape anisotropy favors in-plane magnetization direction.

We have next performed a more detailed local analysis of MAE for a smaller Co nanocrystal made of 110 atoms (shown on the right panel of Fig. 2). For such a relatively small crystal, ab initio DFT calculations within FT approach can be also carried out and compared with TB results. Fig. 7 reports atom-resolved MAE for such pyramid. The atoms of each atomic layer are numbered starting from the corner and going anticlockwise along the spiral to the centre of the plane, as shown in Fig. 6 (a) for the base layer. The other layers are numbered in the same way. Again, a qualitatively good agreement has been found between TB and DFT calculations. Interestingly, we found a sign change of MAE between atomic layers: the MAE favors in-plane magnetization for the first and third layers and out-of-plane magnetization for the mid- To get more insight into the local composition of MAE, we have looked at its distribution in the real space as defined in Eq. 5. Such a real space representation of MAE for the previously studied 110 atoms Co pyramid is shown in Fig. 7. Interestingly, there are regions of both positive as well as negative MAE around each atom (Fig. 7 a, b), in relative proportion which changes from layer to layer. This leads, on average, to the change of sign for atomic MAE vs. the layer observed in Fig. 6. We notice moreover that positive and negative regions of MAE have different spatial localization: while the first one extends out of atomic planes (along the z axis) the second one is mostly localized in the xy plane. This observations can be important when studying the MAE modification due to deposition of pyramids on various substrates.

IV. Conclusion

We have presented a combined TB and DFT study of magnetocrystalline anisotropy of iron (bcc) and cobalt (fcc) slabs and nanocrystals. The nanocrystals are in shape of truncated pyramids with the same length to height ratio as in the experiment. Thanks to the use of the Force Theorem that we have recently implemented in the QE package we have been able to perform a careful local analysis of the MAE in these nanostructures. The TB model is in good agreement with the DFT calculations and gives us confidence in the validity of our TB results for large crystals that cannot be done within the DFT approach. We found a large in-plane anisotropy for Co(001) and a relatively small out-of-plane for Co(111) due to cancelation from the sub-surface layer in the latter case. This is in contrast with iron surfaces since Fe(001) shows a clear out of plane anisotropy. The densest surface shows however a rather small anisotropy for both elements.

These results could have a direct consequence on the magnetic stability of Fe and Co nanocrystals. Indeed, the total MAE is of the same order of magnitude for both Fe and Co nanocrystals, but opposite in sign. This means that while the spin moment of Fe nanocrystals is fixed along the easy outof-plane axis and needs to overcome the high MAE barrier to reverse from positive to negative direction, the magnetic moment of Co nanocrystals is allowed to rotate almost freely (with a very low in-plane anisotropy barrier) in the easy basal plane. One can thus conclude that Fe nanoclusters should be better candidates for magnetic storage applications. Our local analysis, however, indicates that the MAE of nanocrystals could be substantially altered, for instance, by their covering with a mono-layer of another chemical element or by their deposition on various substrates (SrTiO 3 (001), Cu, Au etc).

  FIG. 1: (color online) (a) Co deposition onto a 350 • C SrTiO3(001)c(4×2) substrate followed by a 320 • C anneal gives rise to truncated pyramid shaped nanocrystals as shown in the STM image, 80×80 nm 2 , Vs = +1.0 V, It = 0.1 nA. (b) height to length ratio constant of Co nanocrystals.

FeFIG. 2 :

 2 FIG. 2: (color online) Examples of truncated-pyramid shaped Fe and Co nanocrystals studied in the present work. The crystals are made of bcc Fe and fcc Co with two types of facets: (001) and (110) for Fe and (001) and (111) for Co, respectively. Their possible size and shape is controlled by length-to-height ratio, l/h, kept to ∼ 1.0 (Fe) and 1.41 (Co) which are close to experimental values, ∼ 1.20 (Fe)16 and ∼ 1.48 (Co). The z axis was chosen to be normal to the pyramid base and the spin moment is rotated in the xz plane forming the angle θ with the z axis.

FIG. 3 :

 3 FIG.3:(color online) Total MAE per unit cell of N atoms (in meV), namely, E ⊥ tot -E tot , versus the Co film thickness N for fcc Co slabs. Squares and circles are for (001) and (111) slabs respectively. TB calculations (blue) are compared with ab initio DFT-GGA fully relativistic calculations (red) until N = 10. Lines are guide for the eyes.

FIG. 4 :FIG. 5 :

 45 FIG.4:(color online) Layer-resolved MAE per Co atom (in meV) of slabs with 20 atomic-layer-thick calculated by Tight-Biding (top) and DFT-GGA (bottom) within force theorem approximation. Blue squares and red circles are for (001) and (111) slabs, respectively. Lines are guide for the eyes.

FIG. 6 :FIG. 7 :

 67 FIG. 6: (color online) Atom-resolved MAE for Co nanocrystal made of 110 atoms: (a) trajectory for numbering the base layer atoms starting from the corner and going along the spiral to the center. The atoms of other layers are numbered in the similar way. (b) MAE per atom in meV within TB approach. (c) MAE per atom in meV from DFT-GGA calculations.

TABLE I :

 I TB results: Local analysis of MAE for

Acknowledgement

The research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP7/2007-2013) / ERC grant agreement n • 259297. The authors would like to thank the Royal Society and DSTL for funding and Chris Spencer (JEOL UK) for valuable technical support. This work was performed using HPC resources from GENCI-CINES (Grant Nos. x2013096813).