Radhouène Neji

Jean-François Deux

Gilles Fleury

Mezri Maatouk

Georg Langs

Jean-Philippe Thiran

Guillaume Bassez

Alain Rahmouni

Nikos Paragios

A Kernel-based Approach to Diffusion Tensor and Fiber Clustering in the Human Skeletal Muscle

Keywords: DTI, Diffusion tensor, Fiber, Kernels, Clustering, Human skeletal muscle

In this report, we present a kernel-based approach to the clustering of diffusion tensors in images of the human skeletal muscle. Based on the physical intuition of tensors as a means to represent the uncertainty of the position of water protons in the tissues, we propose a Mercer (i.e. positive definite) kernel over the tensor space where both spatial and diffusion information are taken into account. This kernel highlights implicitly the connectivity along fiber tracts. We show that using this kernel in a kernel-PCA setting compounded with a landmark-Isomap embedding and k-means clustering provides a tractable framework for tensor clustering. We extend this kernel to deal with fiber tracts as input using the multiinstance kernel by considering the fiber as set of tensors centered in the sampled points of the tract. The obtained kernel reflects not only interactions between points along fiber tracts, but also the interactions between diffusion tensors. We give an interpretation of the obtained kernel as a comparison of soft fiber representations and show that it amounts to a generalization of the Gaussian kernel Correlation. As in the tensor case, we use the kernel-PCA setting and k-means for grouping of fiber tracts. This unsupervised method is further extended by way of an atlas-based registration of diffusion-free images, followed by a classification of fibers based on non-linear kernel Support Vector Machines (SVMs) and kernel diffusion. The experimental results on a dataset of diffusion tensor images of the calf muscle of 25 patients (of which 5 affected by myopathies, i.e. neuromuscular diseases) show the potential of our method in segmenting the calf in anatomically relevant regions both at the tensor and fiber level.

Une approche basée sur des noyaux pour le groupement de tenseurs de diffusion et de fibres dans le muscle squelettique

Résumé : Dans ce rapport, nous présentons une approche basée sur des noyaux pour le groupement de fibres dans le muscle squelettique. En se basant sur l'intuition physique derrière les tenseurs de diffusion comme moyen de représenter l'incertitude sur la position des molécules d'eau dans les tissus, nous proposons un noyau de Mercer (défini positif) sur l'espace des tenseurs qui tient compte aussi bien de l'information spatiale que de l'information de diffusion. Le noyau met l'accent implicitement sur la connectivité le long des fibres. Nous montrons que l'utilisation de ce noyau pour une analyse en composantes principales combinée à l'algorithme Isomap fournit un cadre pratique pour regrouper les tenseurs. Nous étendons ce noyau pour pouvoir segmenter des fibres en considérant une fibre comme étant un ensemble de tenseurs centrés sur les points échantillonés le long de sa trajectoire. Le noyau obtenu reflète aussi bien les intéractions spatiales que les intéractions entre tenseurs de diffusion. Nous donnons également une interprétation de ce noyau comme étant une faon de comparer des représentations probabilistes des fibres. Nous étendons la méthode pour effectuer une segmentation supervisée à l'aide d'un atlas en utilisant des SVMs non linéaires combinés avec une diffusion base sur un champ de Markov. Les résultats expérimentaux effectués sur les données de mollet de 25 patients (dont 5 atteints de myopathies) montrent le potentiel de notre méthode pour segmenter le muscle en régions anatomiquement significatives.

Mots-clés : IRM de diffusion, Tenseur de diffusion, Fibre, Noyaux, Groupement, Muscle squelettique

1

Conclusion

RR n 6686

1 Introduction

Context and Motivation

Diffusion Tensor Imaging (DTI) is a modality that has been used extensively in the study of the connectivity within the different anatomical structures in the human brain [START_REF] Le Bihan | Diffusion tensor imaging: Concepts and applications[END_REF]. The acquisition setting allows to compute a field of 3 × 3 symmetric positive definite matrices that model the uncertainty information (covariance) of a Gaussian distribution over the displacements of water protons in the tissues. DTI has attracted much interest because of its potential in unveiling the fiber tracts in brain white matter, since diffusive transport in organized tissues is more important along fiber directions and is rather hindered in the orthogonal directions. However, little work has been done to harness the potential of this modality to explore other regions where water diffusion can carry rich information about the architecture of the fibers and the underlying structure and organization. The human skeletal muscles and more specifically the lower leg are of particular interest because they present an ordered structure of elongated myofibers. Moreover, the perspective of studying the effect of neuromuscular diseases (myopathies) on water diffusion in the human skeletal muscle and providing a tool for early diagnosis and quantitative assessment of muscular weakness and atrophy due to illness based on DTI information is enticing. For the time being, it is unclear how the muscular diseases affect the myofibers, therefore separating the different muscle groups of the skeletal muscles of the lower leg (calf) in regions consistent with anatomical knowledge is a crucial step for a localized quantitative study and comparison of healthy and diseased fiber bundles.

In the following, we review the existing methods for tensor and fiber clustering.

Previous Work

Diffusion Tensor Clustering

The existing diffusion tensor clustering algorithms can be subdivided in three groups. The first class of methods uses a variational approach along with an adequate distance over the manifold of 3 × 3 symmetric positive definite matrices. For instance in [START_REF] Lenglet | DTI segmentation by statistical surface evolution[END_REF], a 3D surface is evolved using an implicit level set representation to segment a region of interest where the spatial gradient is computed using the geodesic distance and the distributions of the tensors in each region are modeled as Gaussian. In [START_REF] Wang | DTI segmentation using an information theoretic tensor dissimilarity measure[END_REF], the Mumford-Shah functional is minimized using a distance between tensors derived from the Burg divergence. A level set technique is also used in [START_REF] Melonakos | Finsler tractography for white matter connectivity analysis of the cingulum bundle[END_REF] to extract the cingulum based on the Finsler metric. Similarly in [START_REF] Jonasson | A level set method for segmentation of the thalamus and its nuclei in DT-MRI[END_REF], several similarity measures are investigated and guide the evolution of coupled level sets.

The second class of methods uses common clustering algorithms to achieve segmentation of the tensor field. In [START_REF] Suyash | A fuzzy, nonparametric segmentation framework for DTI and MRI analysis: With applications to DTI-tract extraction[END_REF], the authors propose to use the Log-Euclidean metric to obtain a kernel density estimate of the probability distribution of tensors and include it in a fuzzy k-means framework where the spatial interactions are handled using local Gaussian kernels INRIA with a fixed bandwidth. In [START_REF] Duan | Thalamus segmentation from diffusion tensor magnetic resonance imaging[END_REF], mean-shift clustering is applied for the segmentation of the Thalamus using Gaussian kernels both in the tensor and in the position space.

The third class of methods consists of graph theoretical approaches and manifold learning techniques that try to capture the structure of the tensor field and propagate the information between neighbors. In [START_REF] Yonas | DT-MRI segmentation using graph cuts[END_REF], a graph-cut approach is used with seed-point initialization. Spectral clustering is performed in [START_REF] Ziyan | Segmentation of thalamic nuclei from DTI using spectral clustering[END_REF] through the eigenanalysis of an affinity matrix between tensors based on a selected similarity measure. In [START_REF] Goh | Clustering and dimensionality reduction on Riemannian manifolds[END_REF], several manifold learning techniques are extended to deal with non-Euclidean spaces and spatial connectivity is ensured by considering isotropic neighborhoods. The particular case of Locally Linear Embedding (LLE) is further discussed in [START_REF] Goh | Segmenting fiber bundles in diffusion tensor images[END_REF] with different choices of tensor metrics.

Fiber Clustering

Graph theoretical approaches have been particularly popular for fiber clustering. For instance, a Gaussian model of the fibers and a normalized-cut approach based on the Euclidean distance between the moments is presented in [START_REF] Brun | Clustering fiber tracts using normalized cuts[END_REF]. In [START_REF] Odonnell | White matter tract clustering and correspondence in populations[END_REF], spectral clustering along with the Hausdorff distance between fibers is considered. The method presented in [START_REF] Brun | Coloring of DT-MRI fiber traces using Laplacian eigenmaps[END_REF] relies on Laplacian Eigenmaps and similarity between fibers is determined using their end points. In [START_REF] Tsai | Fiber tract clustering on manifolds with dual rooted-graphs[END_REF], the authors suggest another manifold learning technique by constructing a graph-based distance that captures local and global dissimilarities between fibers and use LLE for clustering of the tracts. Curve modeling has attracted attention and was handled in [START_REF] Maddah | A unified framework for clustering and quantitative analysis of white matter fiber tracts[END_REF] by defining a spatial similarity measure between curves and using the Expectation-Maximization algorithm for clustering. More recently, fibers were represented in [START_REF] Savadjiev | Streamline flows for white matter fibre pathway segmentation in diffusion MRI[END_REF] using their differential geometry and frame transportation and a consistency measure was used for clustering Of particular interest in the field are the supervised methods that try to achieve a segmentation consistent with a predefined atlas. Registration of B0 images and a hierarchical classification of fibers is performed in [START_REF] Maddah | Automated atlas-based clustering of white matter fiber tracts from DTMRI[END_REF] using the B-spline representation of fibers. The method proposed in [START_REF] Odonnell | White matter tract clustering and correspondence in populations[END_REF] is further extended in [START_REF] Odonnell | Automatic tractography segmentation using a high-dimensional white matter atlas[END_REF] by means of a Nystrom approximation of the out-of-sample extension of the spectral embedding.

Contributions of this Work

We address several issues in this report. While the previous fiber clustering methods discard the tensor information and rely on the obtained tracts, the method we propose handles tensors and is easily extended to segment fibers while taking into account the initial tensor field. We also bridge the gap between tensor and fiber clustering. To this end, we propose a kernel over the tensor space that is consistent with the physical intuition of diffusion tensors as representing the covariance of the probability distribution of water protons positions. Unlike the existing similarity measures and the use of isotropic Gaussian kernels or neighborhoods for spatial interaction, the proposed kernel quantifies not only the dissimilarity between tensors, but takes into account their localization in space in a tractable way and enhances implicitly the connectivity along fiber tracts, i.e. in the feature space provided by the kernel embedding, tensors which are aligned will be closer than tensors which do not lie on the same fiber tract, which is not guaranteed by spatially isotropic Gaussian kernels. We use kernel-PCA for embedding in a Euclidean space and landmark-Isomap for information propagation along diffusive pathways to obtain the final embedding. This is important as it allows to reflect the diffusion flow along the muscle fibers which are characterized by an elongated structure. The clustering is done in the embedding space by a plain k-means algorithm.

The proposed kernel is extended to deal with fiber tracts as input of the clustering algorithm by way of the summation kernel, which is a handy way to define kernels over sets. Not only we take into account in this kernel the interactions between the points (as spatial positions) but also the information provided by the whole tensor field. Given that the fiber pathways are provided, only kernel-PCA and k-means are sufficient for clustering in the fiber domain. We give an interpretation of the fiber kernel as a comparison of soft representations of the fiber tracts and show that it provides a natural generalization for Gaussian kernel Correlation.

In a last step, we show how this unsupervised algorithm can be made supervised in an efficient way. Unlike the work in [START_REF] Odonnell | Automatic tractography segmentation using a high-dimensional white matter atlas[END_REF], given that the kernel satisfies the Mercer conditions, we can look at the problem from a discriminative perspective and avoid the Nystrom outof-sample extension by using a kernel SVM classification and Markov Random Field kernel diffusion of the obtained classification scores in order to segment the calf muscle in regions consistent with a previously defined atlas.

The remainder of this report is organized as follows: In section 2, we discuss the kernel over the tensor space and provide an analysis of its properties and advantages. We also discuss its use in a kernel-PCA and landmark-Isomap setting for clustering of diffusion tensors. Based on the summation kernel, we provide an extension of this kernel to deal with fiber tracts in section 3, where we also discuss the interpretation of the obtained kernel as a means to measure similarity between soft representations of the fiber tracts. We discuss also the special case of Gaussian kernel Correlation. In section 4, we provide the supervised version of the algorithm based on kernel SVMs and kernel diffusion. Section 5 is dedicated to the experimental results and we discuss the perspectives of this work in section 6.

A Probability Kernel on Tensors

Diffusion tensors measure the motion distribution of water molecules. More explicitly, they refer to the covariance of a Gaussian probability over the displacements r of the water protons given a diffusion (mixing) time t:

p(r|t, D) = 1 det(D)(4πt) 3 exp - r t D -1 r 4t (1)

INRIA

Given a diffusion tensor D localized at voxel x, we can obtain the probability of the position y of the water molecule previously localized at x in a straightforward way:

p(y|x, t, D) = 1 det(D)(4πt) 3 exp - (y -x) t D -1 (y -x) 4t (2)
Therefore, a natural way to define a kernel over the tensor space where position is taken into account is to consider the expected likelihood kernel [START_REF] Jebara | Probability product kernels[END_REF]. Let us consider two tensors D 1 and D 2 localized at x 1 and x 2 respectively, and a diffusion time t. The expected likelihood kernel

k t ((D 1 , x 1 ); (D 2 , x 2
)) between the pairs (D 1 , x 1 ) and (D 2 , x 2 ) is defined as the expectation of Gaussian probability p 2 (y|x 1 , t, D 1 ) under the probability law of p 1 (y|x 2 , t, D 2 ) and is given by the following expression:

k t ((D 1 , x 1 ); (D 2 , x 2 )) = E p2(y|x2,t,D2) (p 1 (y|x 1 , t, D 1 )) (3) = p 1 (y|x 1 , t, D 1 )p 2 (y|x 2 , t, D 2 )dy (4)
Note that the diffusion time t is naturally a parameter for this kernel. Using the expression provided in 2, one can obtain the following closed-form expression of this kernel:

k t ((D 1 , x 1 ); (D 2 , x 2 )) = 1 (4πt) 3 k 1 (D 1 , D 2 )k 2 ((D 1 , x 1 ); (D 2 , x 2 )) (5)
where

k 1 (D 1 , D 2 ) = 1 det(D 1 + D 2 ) k 2 ((D 1 , x 1 ); (D 2 , x 2 )) = exp - 1 4t (x t 1 D -1 1 x 1 + x t 2 D -1 2 x 2 ) × exp 1 4t (D -1 1 x 1 + D -1 2 x 2 ) t (D -1 1 + D -1 2 ) -1 (D -1 1 x 1 + D -1 2 x 2 ) (6)

Properties of the Tensor Kernel

The kernel stems from an L 2 inner product defined on the Hilbert space of square-integrable functions, to which Gaussian probability densities belong. Therefore the kernel verifies the Mercer conditions, i.e. it is positive definite. We hereafter provide an analysis of this kernel:

The first term k 1 (D 1 , D 2 ) may be rewritten as follows:

k 1 (D 1 , D 2 ) = 1 det(D 2 ) 1 det(Id + D 1 D -1 2 ) = 1 det(D 1 ) 1 3 i=1 (1 + λ i ) (7)
RR n 6686

where Id is the 3 × 3 identity matrix and λ i are the generalized eigenvalues of the pair of matrices (D 1 , D 2 ). This is reminiscent of the geodesic distance on the manifold of 3 × 3 symmetric positive definite matrices d = 3 i=1 (log(λ i )) 2 [START_REF] Pennec | A Riemannian framework for tensor computing[END_REF] which is also based on the generalized eigenvalues: the distance (respectively the kernel) is increasing (respectively decreasing) with increasing (respectively decreasing) generalized eigenvalues, which is a reasonable behavior (recall that the kernel reflects similarity). The symmetry is ensured in the geodesic distance by the squared logarithm function because the latter is invariant with respect to the inverse transformation λ i → 1 λi (recall that the generalized eigenvalues of the pair (D 2 , D 1 ) are the inverse of those of the pair (D 1 , D 2 )), while it can be seen that the factor

1 √ det(D2)
has a similar role since it also preserves the symmetry. Note that the original expression is clearly symmetric.

Two special cases of the second factor k 2 are interesting:

1. When D 1 = D 2 = D k 2 ((D, x 1 ); (D, x 2 )) = exp - 1 8t (x 1 -x 2 ) t D -1 (x 1 -x 2 ) (8)
As expected, when the tensors are equal, what appears is the Mahalanobis distance between positions x 1 and x 2 with respect to D. In particular when the tensor

D is isotropic, i.e. D = µId, k 2 ((D, x 1 ); (D, x 2 )) = exp -1 8µt ||x 1 -x 2 )|| 2
, which is plainly a Gaussian kernel between the points x 1 and x 2 . 2. When

x 1 = x 2 = x, k 2 ((D 1 , x); (D 2 , x)) = 1,
which means that the kernel k reduces to k 1 . Again this was expected since the kernel will rely only on tensor similarity if there is no difference in spatial positions.

The first special case is of particular interest in diffusion tensor analysis to enhance the connectivity between tensors which are aligned on the same fiber tract. Let us consider the tensor configuration in figure 1 where all tensors are equal to D = µ( The tensors are all equal yet the second term k 2 allows to affect more affinity between tensors 1 and 2 than between tensors 1 and 3. Indeed, we can compute the kernel values to obtain

k 2 ((D, x 1 ); (D, x 2 )) = exp - d 2 8νt ( 9
)
k 2 ((D, x 1 ); (D, x 3 )) = exp - d 2 8µt ( 10
)
where

d = ||x 1 -x 2 || = ||x 1 -x 3 ||. Since ν > µ,
we can see that the similarity between the tensors 1 and 2 is higher than the similarity the tensors 1 and 3, despite the fact that the tensors are all equal.The kernel captures locally the fiber structure.

INRIA

Embedding of the Tensors through Kernel PCA and Landmark Isomap

Let us consider the N pairs (D i , x i ) i=1...N representing a tensor field. We construct the

N × N kernel matrix K of entries K ij = k t ((D i , x i ); (D j , x j
)) for a fixed diffusion time t and normalize it to obtain K such that

K ij = Kij √ KiiKjj
. These pairs are then embedded in a k-dimensional Euclidean space using an eigenvalue decomposition of K = USU t where U is an orthogonal N × k matrix and S is a k × k diagonal matrix. The coordinates of the embedded tensors are given by the N × k matrix X = U √ S where √ S is obtained by setting the diagonal elements of S to their square roots [START_REF] Scholkopf | Nonlinear component analysis as a kernel eigenvalue problem[END_REF]. Each row m of X holds the coordinates in the feature space of the m-th pair.

Given the k-dimensional representation X of the tensor field, one has to propagate the local interaction between the tensors and take into account the distance along diffusive pathways, i.e. simulate the water flow along these trajectories. This is done using the Isomap algorithm which is based on two steps [START_REF] Tenenbaum | A global geometric framework for nonlinear dimensionality reduction[END_REF]:

1. Between every two points of the dataset, find the shortest path using the Dijkstra algorithm based on the Euclidean distance in the feature space and compute a new distance matrix that holds the lengths of the path. This step is important in diffusion tensor analysis in order to propagate diffusion information along fiber tracts.

2. Perform Multidimensional Scaling (MDS) to obtain the new embedding, i.e. a configuration of points that respects approximately the distance matrix computed in the previous step. Since the kernel enhances fiber connectivity, we expect the new configuration to reflect the diffusion flow in the tissues.

Neji & al.

Note that in practice we use a faster version of this algorithm called landmark-Isomap [START_REF] De | Global versus local methods in nonlinear dimensionality reduction[END_REF] that reduces the computational time of the first step by computing the distance of the points to a reduced set of landmarks chosen randomly in the dataset. The clustering is done afterwards using a plain k-means algorithm. Note that the kernel PCA step amounts to denoising in the feature space and that from a theoretical point of view (with perfect data), we could have used from the outset the Euclidean distance (L 2 norm) implied by the kernel k t in the Isomap algorithm without the kernel PCA projection, since the latter requires only the distances between points given by the following expression:

d = k t ((D 1 , x 1 ); (D 1 , x 1 )) + k t ((D 2 , x 2 ); (D 2 , x 2 )) -2k t ((D 1 , x 1 ); (D 2 , x 2 )) (11)
In the following section, we extend in a tractable way the proposed kernel defined over the tensor field to the fiber domain.

Probability Kernel in the Fiber Domain

The fiber trajectories are obtained through the integration of the vector field of principal directions of diffusion. Based on the continuous tensor field approximation (by means of interpolation), we represent each fiber tract as a sequence of tensors localized in spatial positions, i.e. is a set of pairs τ i = (D i , x i ) i=1...n where n is the number of points lying on the fiber. Note that the tractography already requires tensor interpolation and that the interpolated tensors are therefore kept for kernel computation. So it is natural to extend the tensor kernel using a kernel over sets. We simply use the summation kernel [START_REF] Haussler | Convolution kernels on discrete structures[END_REF] to obtain the following Mercer kernel K t between two fibers F 1 and F 2 :

K t (F 1 , F 2 ) = 1 n 1 1 n 2 (Di,xi)∈F1 (Dj ,xj )∈F2 k t ((D i , x i ); (D j , x j )) (12)
where n 1 (resp. n 2 ) is the number of points of the fiber F 1 (resp. F 2 ). This kernel sums the interactions between tensors belonging to the fiber tracts. It captures the diffusion and spatial links between diffusive pathways. It is important to notice that while all the interactions are summed, the diffusion time t acts as a scale parameter. Therefore for a suitable choice of t, a tensor interacts only with tensors lying in a local neighborhood and far-away tensors have a negligible impact on the summation. As in the case of tensors, the segmentation of the fiber tracts is achieved using kernel PCA and k-means clustering.

A Physical Interpretation of the Fiber Kernel

One can see that the summation kernel is simply the expected likelihood kernel between distributions providing a soft representation for fibers. More explicitly, we consider a dynamical system where a particle can be initially at a position x i on the fiber tract and moves

INRIA

to a position y with the following probability

p(y|t, (D i ) i=1...n ) = n i=1 p(x i )p(y|x i , t, D i ) (13)
With a uniform prior distribution on x i , p(y|t,

(D i ) i=1...n ) = 1 n n i=1 p(y|x i , t, D i ).
If the initial positions were independent (which is not the case because they are the result of the tractography), this would have amounted exactly to an adaptive kernel density estimation of the position of the water molecules along the fiber tracts where the point-dependent Gaussian kernels use the diffusion tensors as covariance matrices to model the uncertainty. However the distributions in equation 13 still provide a soft representation of the fibers and measure the compactness of the spatial configuration of the fiber tract. By bilinearity of the expected likelihood kernel, it is straightforward to see that the expected likelihood kernel of the distributions given in Equation 13is exactly the summation kernel.

Links with Gaussian Kernel Correlation

When considering the special case where the tensor field is constant (equal to a tensor D), we obtain

K t (F 1 , F 2 ) ∝ 1 n 1 1 n 2 xi∈F1 xj ∈F2 exp - 1 8t (x i -x j ) t D -1 (x i -x j ) (14)
which is the anisotropic form of the Gaussian kernel correlation [START_REF] Tsin | Kernel Correlation as an Affinity Measure in Point-Sampled Vision Problems[END_REF]:

KG t (F 1 , F 2 ) ∝ 1 n 1 1 n 2 xi∈F1 xj ∈F2 exp - ||x i -x j )|| 2 8µt ( 15
)
where D = µId which means that we suppose the tensors are isotropic. We can see that the proposed kernel deals with a generic tensor field and provides a generalization of the Gaussian kernel correlation. A by-product of this reasoning is that Gaussian kernel correlation is a kernel on fibers that considers only point positions. One could have seen from the outset that is a Mercer kernel since it is a summation of Mercer (Gaussian) kernels. This will be particularly useful to learn spatial interactions between fibers as it will be detailed in the next section.

A Supervised Fiber Clustering Framework 4.1 Supervised Kernel SVM Learning of Atlas Fibers

Given an atlas of fibers segmented by an expert in R regions, we can learn the spatial interactions between these fibers using the Gaussian kernel correlation in Equation 15. Indeed, it can be used as an input in a kernel Support Vector Machines (SVMs) [START_REF] Vapnik | Statistical Learning Theory[END_REF] to learn boundaries between the different segmented regions, including background fibers that are not localized in the region of interest. The kernel SVMs provide support vectors (fibers), which are the fibers that define the decision boundaries. Note that the SVMs are used in a one-against-one fashion in order to deal with multiple regions. A point of interest here is that from the initial set of training fibers, we have only to keep a sparse subset of support fibers that will guide the classification process. The fact that the Gaussian kernel correlation defines a kernel allows to avoid the Nystrom approximation of the out-of-sample extension of the spectral embedding proposed in [START_REF] Odonnell | Automatic tractography segmentation using a high-dimensional white matter atlas[END_REF] by looking to the problem from a discriminative perspective as opposed to a projection approach.

Given a new set of fibers that we wish to segment in a manner consistent with the already-defined atlas, we start by finding the affine transformation that maps the diffusionfree (B0) image of the training atlas to the corresponding one in the testing dataset, as in [START_REF] Maddah | Automated atlas-based clustering of white matter fiber tracts from DTMRI[END_REF]. This tranformation is subsequently used to register the testing set of fibers to the space of the support fibers obtained from the atlas. We are therefore able to compute the scores of the R(R-1) 2 pairwise SVM classifiers on the testing dataset. In the following subsection, we show how to use the SVM classification results to obtain the segmentation while respecting the interactions provided by the fiber kernel.

Markov Random Field Segmentation of the Fiber Tracts

We start by embedding the fibers that are to be segmented in a k-dimensional space using kernel PCA. Then for each fiber, we consider its k-nearest neighbors according to the Euclidean distance in the feature space. Based on the defined spatial neighborhoods, we use a Markov random field model and minimize the following energy:

E = Fi u s (l(F i )) + λ Fi,Fj ∈N (Fi) u p (l(F i ), l(F j )) (16)
where l(F i ) is the classification label of the fiber F i and N (F i ) its neighborhood. The first term is provided by the classifiers (data term) and set to u s (l

(F i )) = -log( 2k l R(R-1) )
where k l is the number of pairwise classifiers voting for class l (if k l = 0 , we choose u s (l(F i )) =log(ǫ), where ǫ is a relatively low value).

The second term is used for kernel-based regularization and set to u p (l

(F i ), l(F j )) = K t (F i , F j )(1 -δ(l(F i ), l(F j )))
where δ is the Kronecker delta. The constant λ is a trade-off parameter between the two terms. This amounts to a kernel graph diffusion of the initial labels provided by the SVM classifiers and accounts for the regularity of the segmentation. The energy is optimized using the algorithm described in [START_REF] Komodakis | Fast, approximately optimal solutions for single and dynamic MRFs[END_REF].

Results and Experiments

The human skeletal muscles and more specifically the lower leg are of particular interest because they present an ordered structure of elongated myofibers where some muscle groups INRIA differ only subtly in their direction. In order to study the effect of neuromuscular diseases (myopathies) on water diffusion in the muscles, segmentation in regions consistent with anatomical knowledge is a crucial preliminary step before a localized quantitative study of diffusion properties in the fiber bundles for healthy and diseased tissues. Diffusion-based studies of the human skeletal muscle [START_REF] Galban | Diffusive sensitivity to muscle architecture: a magnetic resonance diffusion tensor imaging study of the human calf[END_REF][START_REF] Sinha | In vivo diffusion tensor imaging of human calf muscle[END_REF][START_REF] Damon | Validation of diffusion tensor mri-based muscle fiber tracking[END_REF] focused on studying variation across subjects of scalar information like diffusivity (trace), fractional anisotropy, pennation angles (the orientation of muscular fibers with respect to tendons), etc. A supervised SVM tensor classification algorithm was proposed in [START_REF] Neji | Support vector driven Markov random fields towards DTI segmentation of the human skeletal muscle[END_REF].

Twenty-five subjects (twenty healthy patients and five patients affected by myopathies) underwent a diffusion tensor imaging of the calf muscle using a 1.5 T MRI scanner. The following parameters were used : repetition time (TR)= 3600 ms, echo time(TE) = 70 ms, slice thickness = 7 mm and b value of 700 s.mm -2 with 12 gradient directions and 13 repetitions. The size of the obtained volumes is 64 × 64 × 20 voxels with a voxel resolution of 3.125 mm × 3.125 mm × 7 mm. We acquired simultaneously high resolution T1-weighted images that were segmented manually by an expert into 7 muscle groups to provide the ground truth and fiber trajectories were reconstructed using [START_REF] Fillard | Medinria: DT-MRI processing and visualization software[END_REF]. To give an idea about the muscle architecture in the calf, we present in [Fig. 2 (a)] a manual segmentation overlaid on an axial slice of a high-resolution T1-weighted image. The following muscle groups are considered: the soleus (SOL), lateral gastrocnemius (LG), medial gastrocnemius (MG), posterior tibialis (PT), anterior tibialis (AT), extensor digitorum longus (EDL), and the peroneus longus (PL).

In the following , we present the obtained experimental results on a synthetic dataset and for tensor classification and both supervised and unsupervised fiber bundling of the lower leg muscles.

Tensor Clustering

Preliminary Experiments

We first generated a 20 × 40 lattice of synthetic tensors composed of two close fiber bundles. The first bundle has a vertical principal direction, the second starts with a vertical direction then deviates with a 45 • angle [Fig. 3 (a)]. We added a Gaussian noise of standard deviation 10 • to these directions. The eigenvalues of the tensors were set to {2 10 -3 , 1.5 10 -3 , 10 -3 }. We tested the following values of diffusion time t: {10 4 , 10 5 }. We compare the behavior of the kernel PCA + Isomap embedding with spectral clustering using the metric

d(D 1 , D 2 ) = arccos(< → e 1 , → e 2 >)
where → e 1 (resp. → e 2 ) is the principal direction of diffusion of D 1 (resp. D 2 ), as in [START_REF] Ziyan | Segmentation of thalamic nuclei from DTI using spectral clustering[END_REF]. Following [START_REF] Ziyan | Segmentation of thalamic nuclei from DTI using spectral clustering[END_REF], the scale parameter in the affinity matrix is set as the sample variance of d between neighboring tensors. The clustering is obtained using k-means with 50 restarts in the spectral embedding space. [Fig. 3 (c)] shows the segmentation result obtained by our approach (stable across the tested set of diffusion times) and [Fig. 3 (b)] shows its counterparts for spectral clustering. We can notice that unlike spectral clustering, the proposed algorithm finds a clustering solution which is more compatible with the tensor

Classification of the Muscle Groups

For each subject, a region of interest (ROI) was manually delineated and we tested the performance of the tensor clustering algorithm both for healthy and diseased subjects. In regions affected by myopathies, the tensors have a relatively small volume since fat replaces the fibers (as can be seen in [Fig. 2 (b)]) and were eliminated through simple thresholding over the determinant. In all the experiments, we set the diffusion time t to t = 2 10 5 and we used a ten-dimensional embedding. For a quantitative evaluation of the method, we follow the validation protocol proposed in [START_REF] Ziyan | Segmentation of thalamic nuclei from DTI using spectral clustering[END_REF]: for all the 25 subjects, the manually delineated ROI was segmented at two different levels (in 7 and 10 classes). The resulting clusters are then classified according to the labels given by the expert. As in [START_REF] Ziyan | Segmentation of thalamic nuclei from DTI using spectral clustering[END_REF], several clusters are allowed to have the same label. We also test the algorithm on a section near the knee which is characterized by a higher amount of noise and artifact obtained automatically by a threshold (set to 20) on the diffusion-free images. We report in [Fig. 6] the boxplots of the dice volume overlap with the expert labels for the 25 patients and its counterpart for the spectral clustering as described in the previous subsection. We can see that our algorithm performs slightly better for the case of the manually-delineated ROI, however for the noisy automatic ROI, spectral clustering is misled by isolated points (see [Fig. 5 (c)]), whereas the performance of our algorithm is not significantly worsened. Note also that the thresholding over the determinant removes some of the tensors originating from the artifact. From a qualitative point of view, one can see in [Fig. 5 (a), (b)] that the algorithm was able to segment correctly fine structures like AT and PL and the segmentation result is rather smooth.

Fiber Clustering

Unsupervised Fiber Clustering

To test the unsupervised kernel-PCA clustering algorithm, we only kept the fibers which have a majority of points lying in the manually delineated ROI. The number of fibers in the different datasets ranged approximately from 1000 to 2500, which makes the computation and eigenanalysis of the kernel matrix achievable in a rather reasonable time. The diffusion time t was set to t = 2 10 4 and we used a ten-dimensional embedding for kernel-PCA.

Figure 7 shows the clustering results in 10 (resp. 7) classes for the fiber tracts of a healthy (resp. diseased) subject. The fibers in [Fig. 7 (a)] correspond to the tensors segmented in [Fig. 5 (a)]. It is interesting to notice that despite the fact that the tractography algorithm was unable to recover fiber tracts in the diseased regions due to the presence of degenerate tensors, the clustering algorithm could still segment the fiber tracts of the healthy region in anatomically relevant subgroups. For quantitative assessment, we report in [Fig. 8 (a)] the boxplot of the dice overlap measures of the fiber segmentation with the expert labeling for the 25 subjects for 7 and 10 classes. Overall, the algorithm performs well in separating the regions of the calf muscle with a mean dice coefficient of 79.5% (respectively 80.93%) and a standard deviation of 5.04% (respectively 5.14%) for 7 (respectively 10) classes.

Supervised Fiber Clustering

A manually segmented volume was used as an atlas. Atlas fibers were assigned to a class based on a simple voting procedure: each fiber is classified according to the majority vote class of the voxels it crosses. We experimented both with affine [START_REF] Fillard | Medinria: DT-MRI processing and visualization software[END_REF] and deformable [START_REF] Glocker | Dense image registration through mrfs and efficient linear programming[END_REF] (a) (b) registration to map the B0 images of the testing case to the B0 images of the atlas. We use kernel SVM classification to learn the fibers of the atlas as explained in section 4.2, using 21 one-against-one pairwise classifiers. The scale parameter in the Gaussian correlation kernel was set to σ = 2 √ µt = 10. We report in [Fig. 8 (b)] the boxplot of the dice overlap coefficients both for deformable and affine registration. We can note that we obtain significantly better results with deformable registration, which was expected given the relatively high interpatient variability and that muscles are soft tissues, so the anatomy and shape are likely to vary significantly across patients. In [Fig. 9 (a)], we show an example of supervised segmentation compared with the ground truth in [Fig. 9 (b)]. We can observe that as opposed to the unsupervised setting ([Fig. 7 (a)]), the MG is not oversegmented.

Conclusion

In this report, we proposed a kernel-based method for clustering of both tensors and fibers in diffusion tensor images. It exploits the physical interpretation behind the modality and offers a unified approach towards tensor and fiber grouping. The kernel defined over the tensor space encompasses both localization and diffusion information and naturally reflects tensor alignment along fiber tracts. We showed its flexibility by extending it to deal with fibers and gave the physical intuition behind its mathematical definition as a kernel over sets of tensors. We also showed how to include expert knowledge by means of kernel SVMs.

Future research will focus on the use of the defined kernels for (affine or deformable) registration of diffusion tensors and will explore the possible extensions of the fiber kernel using stochastic process theory, in particular Gaussian processes. The motivation of the work is to use diffusion tensor imaging in order to detect and monitor the progression of skeletal muscle diseases like myopathies.

4 A

 4 Supervised Fiber Clustering Framework 4.1 Supervised Kernel SVM Learning of Atlas Fibers . . . . . . . . . . . . . . . . 4.2 Markov Random Field Segmentation of the Fiber Tracts . . . . . . . . . . . . 5 Results and Experiments 5.1 Tensor Clustering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1.1 Preliminary Experiments . . . . . . . . . . . . . . . . . . . . . . . . . 5.1.2 Classification of the Muscle Groups . . . . . . . . . . . . . . . . . . . . 5.2 Fiber Clustering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2.1 Unsupervised Fiber Clustering . . . . . . . . . . . . . . . . . . . . . . 5.2.2 Supervised Fiber Clustering . . . . . . . . . . . . . . . . . . . . . . . .

  i=1...3 are the canonical basis of R 3 and ν > µ are the eigenvalues of D. The tensors have therefore a principal direction of diffusion along → e 3 .

Figure 1 :

 1 Figure 1: A configuration where the tensor kernel implicitly puts more weight on the connection between D 1 and D 2 than between D 1 and D 3 , reflecting their alignment.

Figure 2 :Figure 3 :Figure 4 :

 234 Figure 2: (a) An example of a manual segmentation of an axial slice of a high-resolution T1-weighted image showing different muscle groups in the calf. (b) An axial slice of a T1weighted image of the calf of a diseased patient where the zone in hypertension is fat that replaced the muscle.

Figure 5 :Figure 6 :

 56 Figure 5: Axial slices of the segmentation of the tensors for (a) a healthy subject in 10 classes, manual ROI (b) a diseased subject in 7 classes where the MG is partially affected, manual ROI (c) noisy automatic ROI of a section near the knee using spectral clustering (d) noisy automatic ROI using our method

Figure 7 :Figure 8 :

 78 Figure 7: Axial, coronal and saggital views of fiber segmentation for (a) a healthy subject in 10 classes, (b) a diseased subject in 7 classes.

Figure 9 :

 9 Figure 9: Axial, coronal and saggital views for (a) supervised classification in 7 classes (b) the ground truth segmentation.

  Introduction 1.1 Context and Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Previous Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.1 Diffusion Tensor Clustering . . . . . . . . . . . . . . . . . . . . . . . . 1.2.2 Fiber Clustering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Contributions of this Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2 A Probability Kernel on Tensors 2.1 Properties of the Tensor Kernel . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Embedding of the Tensors through Kernel PCA and Landmark Isomap . . .

3 Probability Kernel in the Fiber Domain 3.1 A Physical Interpretation of the Fiber Kernel . . . . . . . . . . . . . . . . . . 3.2 Links with Gaussian Kernel Correlation . . . . . . . . . . . . . . . . . . . . .

RR n 6686


Antoine Zink
email: antoine.zink@culture.gouv.fr

Elisa Porto

Methods and Applications of Absolute Chronology LUMINESCENCE DATING OF THE TANAGRA TERRACOTTAS OF THE LOUVRE COLLECTIONS

Keywords: LUMINESCENCE DATING, TANAGRA, DOSE RECONSTRUCTION

 Montreal, 2004), we investigated the authenticity of some 140 terracottas from the Louvre collections by luminescence dating (thermoluminescence and optically stimulated luminescence). For these terracottas being pieces from a collection, only the dose due to internal radioactivity could be derived from the direct measurement of the internal concentration of radioelements. To estimate the contribution of the external environment, we had to reconstruct the moisture content and the gamma dose rate. The moisture content was based on some indications about the site found in 19th century literature. We reconstructed the external dose rate received by the Tanagras from the ratio of the measured palaeodose to 2300 years. This result agrees well with a low radioactive soil known to be characteristic of the mainland Greece. Applying such values of the external dose rate to our Tanagra, we obtained an age distribution showing two maximums, one at 2300 ± 400 years, the other at 150 ± 50 years.

INTRODUCTION

At the end of the 1860.s, the inhabitants of Vratsi in Beoatia (Greece) working in their fields excavated ancient tombs corresponding to the necropolis of the antique city of Tanagra. The ages of the tombs covered a large chronological range from the Geometric (9th century BC) to the Hellenistic (3rd-1st century BC) period [START_REF] Rayet | Les figurines de Tanagra au Louvre (Tanagra figurines of the Louvre)[END_REF]. Small terracotta figurines were found in great numbers inside and outside the tombs of the the Hellenistic period.

The aspect of the Hellenistic figurines met the taste of 19th century middleclass for the reality of the life in the Greek antiquity, as opposed to Gods and Goddesses more appreciated by earlier aristocratic collectors. Within a few years, Tanagra had become a generic term to name a draped young lady. Actually Tanagra figurines not only represent women, a large part of figurines show young boys and actors. In addition, the so-called Tanagras do not come only from Tanagra and its surroundings (Beoatia and Attike). Tanagra-like small figurines were created all over the Mediterranean area (Myrna, Alexandria, Tarent). The Tanagra and Tanagra-like market quickly became a place of speculation [START_REF] Belle | Voyage en Grece (Travel in Grece)[END_REF]. In 1872, the Louvre museum bought one hundred statuettes [START_REF] Rayet | Les figurines de Tanagra au Louvre (Tanagra figurines of the Louvre)[END_REF], followed by other great institutions like the British Museum and by private collectors. This led to the production of fakes. As most of the figurines came from illegal diggings made by peasants far from the eyes of the antiquity merchants, it was easy to produce fakes to cover the large demand [START_REF] Rayet | Les figurines de Tanagra au Louvre (Tanagra figurines of the Louvre)[END_REF]. Three main types of fakes exist. First, the pure inventions of the 19th century, made by contemporary artists, were figurines with a style distinct from the antique figurines, with a multiplication of the folds and draping. Second, were Tanagras obtained by moulding; being proto-industrial products, several generations could be derived from the same prototype. It was then easy to create fakes based on the antique or modern mould. Between the fake and genuine are the pastiches where pieces from various genuine figurines are put together to create a modern piece.

Using thermoluminescence dating, [START_REF] Goedicke | Echtheitsprüfung an Tanagrafiguren nach der Thermolumineszenzmethode (Authenticity of Tanagra figurines based on thermoluminescence methode)[END_REF] showed that some 20% of the Tanagra figurines in the Berlin collection are modern copies. Following the work in Berlin, some twenty-five pieces of the Louvre collection were investigated by thermoluminescence in the Research Laboratory of the French Museums between 1994 and 1997. When in 1998 an exhibition was planned at the Louvre in co-operation with the Montreal Museum of Fine Art, Tanagra -Myth and Archaeology (Paris -Montreal, 2003), a large investigation of some 200 figurines from the Louvre collections was decided. The set included only pieces from Beoatia (Thebes and Tanagra) and Attike (Athens). The aim of the investigation was to characterize the material, clay, and pigment, and to de-fine criterions of authenticity. The present paper reports the result of the luminescence dating work (June 2002[START_REF] Paris | Tanagra: mythe et archéologie (Tanagra, Myth and Archaeology)[END_REF] with the aim to verify the authenticity of the figurines.

EXPERIMENTAL TECHNIQUES Sample preparation

A statuette was sampled using a 1.8 mm diameter tungsten carbide drill. Some 200 mg of powder was collected: 100 mg for luminescence measurements and 100 mg for chemical analysis (ICP-MS / ICP-AES). The choice of the location of the sample on the statuette was made in agreement with the curators. For practical reasons, mainly due to the small size of the pieces, only 140 statuettes were sampled for luminescence dating. We had also to reject the statuettes studied by X-ray radiographies in the past (it was the case of the famous Blue Lady). It should be noted that some sampling was made in the plate supporting the figurine.

After etching in 10% HCl and washing in water, ethanol and acetone, the 4-13 mm fraction was selected by sedimentation into tubes filled with 8 cm height of acetone and deposited on 9.8 mm diameter stainless steel discs. All these operations with the exception of the HCl etching, are automated [START_REF] Zink | Automation of the chemical preparation for luminescence dating[END_REF]. Sixteen discs were prepared for thermoluminescence and four for OSL.

Measurement apparatus

In the 90.s, investigations were done by additive TL using a Daybreak 1100 reader. The luminescence was detected through a combination of 7-59 x BG39 filters by an EMI 9635QA photomultiplier. The heating was limited to 500°C followed by a prompt cooling. The external irradiations were carried out using a Daybreak 801 multi sample irradiator with a 3.7 GBq beta 90Sr/90Y source (3.82 Gy/min on 1/6/2002) and a 18.8 MBq alpha 238Pu source (7.78 mm-2/min on 1/6/2002). For the set of measurements carried out in 2002-2003, the luminescence measurements were performed with a Risø TL/OSL DA-15 set fitted with a 90Sr/90Y source delivering 7.89 Gy/min (on 1/6/2002). All measurements were made using an EMI 9235QA photomultiplier.

We used a double single aliquot regenerated OSL (SAR-OSL) protocol. The irradiations were carried out using the internal 90Sr/90Y source. After each irradiation, the disc was preheated at a temperature T. For each sample we used four discs corresponding to four different temperatures (220°C, 270°C, 320°C and 370°C). The disc was first stimulated by an IR laser diode (830±10 nm; 50 % of 450 mW/cm2 full power) at 60°C during 100s and then stimulated at 125°C for 100s by blue diodes arranged in 21 pairs (470± 30 nm; 50 % of 19 mW/cm2 full power). The luminescence was detected through a 7.5 mm thick U-340 filter. To prevent the leakage of radiation and stimulation to adjacent positions, the discs were placed only on the odd positions. All measurements were recorded by integrating over one second per channel. To determine the alpha efficiency, the disc was irradiated by an alpha external source, and promptly measured by OSL in the same conditions as previously. After normalisation by a test dose, the response was compared to the beta growth curve.

For additive thermoluminescence, the discs were irradiated by the external beta and alpha sources. After irradiation, the discs were stored in darkness at 50°C during 1 to 3 months. After a preheat up to 220°C, the disc was directly cooled to 60°C and then heated up to 650°C and the remaining signal was annealed at 650°C during 2 minutes. The luminescence was detected through a combination of 7-59 x HA-3 Filters. To control fading, a test disc previously annealed at 650°C was irradiated and stored together with the other discs. After measurement, the test disc was irradiated again and promptly (around 4 minutes after irradiation) measured to have a reference.

Alpha counting was carried out with three Daybreak 582/583 alpha counters using a 12 mm diameter cell. Xrays analysis was carried out by energy dispersion spectrometry (EDS) using a scanning electronic microscope.

ICP-AES (Inductively Coupled Plasma -Atomic Emission Spectrometry) and ICP-MS (Mass Spectrometry) were made at the Centre de Recherches Pétrographiques et Géochimiques in Nancy and analysed by A. Bouquillon at the laboratory (Bouquillon et al., in print).

MEASUREMENTS RESULTS Palaeodose

A first test using a SAR-TL protocol was not successful. It was due to a great change in sensitivity after the first measurements (in comparison with the natural signal), which could not be corrected by the test dose response. We attributed this to a modification of the self-absorption of the grains.

The full set of Tanagras was measured by single aliquot regenerated (SAR) OSL and some forty of them were investigated by additive thermoluminescence (TL). No measurable fading effect of the thermoluminescence was observed between one minute and 3 months for our samples. Although we employed a double-SAR OSL protocol, IR-OSL was extremely rare. When we began the study, the a-values were not measured for OSL and we used an average based on the a-value for TL. After a month, we added an easy method to determine this parameter in OSL by adding a cycle in the OSL-SAR protocol with alpha irradiation, measurement and then test dose and measurement. The intensity due to alpha dose was normalised by the test dose. The intensity thus calibrated was plotted on our growth curve. Hence, the beta dose equivalent to the alpha dose was graphically estimated. Later, the avalue was calculated in the same way as in TL. The results had a good precision, better than for TL (Table 1). Both a-values are known only for six samples, and it is difficult to discuss their differences. Obviously the average of the avalues are distinct for TL and OSL. It can correspond to a difference in the physical process. Generally the TL and OSL gave similar equivalent doses. We chose to concentrate our work on the OSL.

Annual dose

The estimation of the annual dose was more complicated. Being collection pieces, only the internal concentration in radioelements could be directly measured giving information on the contributions of alpha and beta dose rate. Nineteen Tanagras had been investigated by alpha counting and EDS in the 90.s. In 2002, all the collection was investigated by ICP-AES (majors) and ICP-MS (traces) to determine the origin of the clays and identify the sites of production. In general, the various methods are in agreement. The only difference between ICPMS and alpha counting comes from the [Th]/[U] ratio. For alpha counting, due to the small size of the sample holder (diameter 12 mm), it is not possible to use pair counting to determine the Thorium content (ca. 1-10 pairs/week). Hence, we used an arbitrary value of 3.167. In fact, the elemental analysis showed a higher [Th]/[U] ratio close to 4.0. Only the concentrations obtained by elemental analysis were used in the age calculation.

Moisture

To estimate the contribution of the external environment, we must reconstruct the moisture content and the gamma dose rate. For practical reasons it was not possible to use a drying oven to measure the moisture inside the Tanagra. Moreover, the water content in the present location of conservation (gallery or storeroom) is not representative of the moisture in the buried conditions. It is necessary to choose an arbitrary value, the most representative to our knowledge about the average moisture percentage. For this purpose we used two indications. Following [START_REF] Aitken | Thermoluminescence Dating[END_REF], the porosity of ceramic varies in the range 5-25 %. The sampling showed that the earth was crumbly. As a result we assumed a high porosity, some 20-25 %. Moreover, the authors from the 19th century described the weakness of the pieces due to high moisture [START_REF] Rayet | Les figurines de Tanagra au Louvre (Tanagra figurines of the Louvre)[END_REF]. Therefore, we suggested a high level of moisture around 16% which corresponds to a fraction of porosity of 60-80 %.

Gamma dose rate

To reconstruct the external dose rate in a general case, we created a database containing more than 725 external dose rates reported in the literature with a worldwide distribution [START_REF] Zink | Bayesian approach applied to authenticity testing by luminescence[END_REF]. More than 95 % of the values are below 2.0 mGy/a. The distribution of the values show two main maximums, the first around 0.4 mGy/a corresponding to low radioactive soils as loess or sediments, the other around 1.2 mGy/a attributed to high radioactive soils such as granites and basaltic rocks. An intermediary maximum around 0.8 mGy/a corresponds to mixed soils, e.g. fluvial deposits from various origins. Continental Greece is mainly constituted by weakly radioactive soils [START_REF] Clouvas | Extended survey of indoor and outdoor terrestrial gamma radiation in Greek urban area by in-situ gamma spectrometry with a portable Ge detector[END_REF]. Hence we expect values close to the first maximum.

Usually, with an artifact of unknown origin (i.e. with an unknown external dose), we assume a local origin of the earth used for the terracotta, the socalled «local hypothesis». Hence, we take for the external gamma dose rate the same value as the internal gamma dose rate. In this case, considering the internal content in radioelements, we found a mean value of 0.62±0.29 mGy/a (Fig. 1). And we assume an uncertainty of 0.6 mGy/a, corresponding to the standard deviation of a uniform distribution between 0 and 2mGy/a. With such uncertainty, higher than the simple distribution of the internal gamma dose rate, all chemical modifications during the manufacture of the piece can be neglected.

For the Tanagras, we used another method. We tried to reconstruct the dose rate experienced by the Tanagras. If the pieces are genuine, their age must be around 2300±100 years. We can assume too, that they laid in only one place, except during the last century. Hence, averaging the palaeodose on 2300 years, we obtain for each Tanagra an annual dose rate. Using the internal dose rate and the fraction of moisture, it is possible to calculate a distribution of the external dose rate among the different statuettes (Fig. 1). This distribution can be fitted by a Gaussian of an average of 0.55 mGy/a and a standard deviation of 0.20 mGy/a.

The two methods yielded similar results for the average external gamma dose rate corresponding to weakly radioactive soils. They gave us the mean ambient dose rate in Beoatia and Attike during the last 2400 years. But for individual Tanagra, the two methods did not always give similar results. When the differences were larger than the uncertainties, we assumed that the local hypothesis failed. The local hypothesis gave us the dose rate of the terracotta, and indirectly of the clays used to produce the piece. On the other hand, the reconstructed dose rate gave us a good estimation of the ambient dose rate, i.e. of the soil where the Tanagras were buried. The agreement between the two hypotheses implied from a radiometric point of view that the object was produced and then conserved in the same region. The rejection of this agreement led us to assume a difference between the place of production and the place of conservation (mainly tombs).

Such analysis could be helpful to identify trade routes between central Greece and volcanic regions like Great Greece or Cyclades. It is known that some Tanagra moulds travelled across the Mediterranean area.

DISCUSSION

For the calculation of the age we decided to take the value for the external dose-rate given by the reconstructed dose rate as 0.55±0.20 mGy/a, which appeared as the most representative of the external dose rate experienced by the Tanagras. We obtained an age distribution showing two maximums, one at 2300±400 years, the other at 150±50 years, corresponding to the age of genuine Tanagras and to the age of the 19th century copies (Fig. 2), respectively. Hence, we were able to separate the two groups. The fakes seem to have undergone similar annual palaeodosimetry as the genuine pieces. In the previous work on the dosimetry in the palaces and museums in Paris and surroundings, we observed a large range of dosimetry from 0.30 to 1.26 mGy/a [START_REF] Castaing | Radiation background due to radioactivity in palaces and museums: influence on TL/ OSL dating[END_REF]. In particular, in the rooms of Greek, Etruscan and Roman Antiquity of the Louvre museum, the observed dose rates span between 0.43 and 1.02 mGy/a. It is difficult to reconstruct the exact dose rate received in the museum as the collections are moved from one room to another.

Our results showed that the same clay groups identified by chemical analysis were observable in the genuine pieces as well as in the fakes, showing that the same clays were used by modern forgers as in antiquity. No modern pigments were used, only the use of Egyptian blue (cuprorivaite) can conversely be of use as a criterion of antiquity. Luminescence dating appears then as the only laboratory method to identify fakes [START_REF] Bouquillon | Authenticité, matières et couleurs. Etude en laboratoire des Tanagréennes du Louvre (Authenticity, maters and colors. Study in laboratory of the Tanagra-like terracotta of Louvre)[END_REF]. After discussion with curators, our results agreed well with their observations. One problem concerned a figurine of an actor (MNB904). A similar figurine, now in the Hermitage Museum (St Petersburg) was actually found in a supervised digging which obviously indicated that this type is ancient. But after discussion, a new examination by the curators showed that the stylistic technique of the piece from the Louvre is in fact fully modern [START_REF] Becq | 236 -Acteur en papposilène parodiant une tragédienne. (Actor as Papposilene parodying a tragic actress)[END_REF].

CONCLUSIONS

To be able to use luminescence dating for pieces from a collection, some hypotheses are needed. For the Tanagras, we have reconstructed the mean ambient dose rate in Beoatia and Attike regions. Our investigation showed that the percentage of fakes in the Louvre collection is small compared with other museums [START_REF] Goedicke | Echtheitsprüfung an Tanagrafiguren nach der Thermolumineszenzmethode (Authenticity of Tanagra figurines based on thermoluminescence methode)[END_REF]. The gap of one year between the purchase of first Tanagras by the Louvre (in 1872) and by Berlin (in 1873) could explain such a difference.

Compared to the simple typo-chronology based on stylistic characteristics, luminescence dating first appears to be less precise. Our results on the Tanagras yield a weighted mean age of 2220±195 years, while curators and art historians can split the Tanagra century (200 . 330 BC) in four or five periods. But the moulding techniques made easy the re-cycling of older models. Some popular figurines could have been produced over more than a century. An ancient mould, or a mould made from ancient pieces, can be used to create modern fakes. In this case, the stylistic definitions based on marble sculptures are not sufficient to ensure authenticity. Then with a lower precision, but a higher reliability, the luminescence dating is safer to date Tanagra figurines.

Fig. 1 .

 1 Fig. 1. Probability density function of the gamma dose rate based on radioelement contents in the figurines (dashed line), compared with the probability of the reconstructed external dose rate (solid line). The reconstructed external dose rate (D.ext) was estimated for each Tanagra using the stylistic age (200 -330 BC), the equivalent dose (De) and the internal dose rate (D.int) (D.ext=De/Age -D.int).

Fig. 2 .

 2 Fig. 2. Probability density function of the age of our Tanagra set (dashed line), based on SAR-OSL measurements, taking 16% for moisture and 0.55±0.20 mGy/ a for external gamma dose rate. The individual age (with 1s uncertainties) are plotted in increased order from bottom to top.

Table 1 .

 1 Alpha efficiency (a-value) measured for OSL and TL. (Lab. No. . laboratory reference number; Inv. No. . Louvre reference number).

	Lab. No. FZ29602	Inv. No. MNB589	A-value OSL 0.210	A-value TL
	FZ29985 FZ30040	MNB595 CA1627	0.160 0.041
	FZ30042 FZ30045	MNB1003 CA3005-2	0.080 0.063	0.080
	FZ30049 FZ30050	MNB481 MNB551	0.110 0.055
	FZ30055 FZ30174	MNB573 MNB585	0.097 0.102	0.120
	FZ30286 FZ30363	MNB489 MNB493	0.061 0.051	0.220
	FZ30364	MNB497
	FZ30365 FZ30366	MNB560 MNB577	0.091 0.060
	FZ30459 FZ30463	MNB1008 MNB492	0.090 0.140
	FZ30464 FZ30467	MNB494 MNB482	0.200 0.120
	FZ30494 FZ30539	MNB559 CA2152	0.080 0.071
	FZ30543 FZ30588	CA2162 CA2168	0.057 0.050
	FZ30589 FZ30590	CA3005-4 CA3005-5	0.059 0.061
	FZ30591 FZ30628	MNB476 MNB596	0.070 0.077	0.130 0.140
	FZ30712 FZ30713	CA588 CA2165	0.020 0.076	0.130
	FZ30750 FZ30751	CA2552 MNB447	0.061 0.160
	FZ30753 FZ30754	MNC730 CA265	0.220 0.140
	FZ31189 FZ31254	CA806 MNC747b	0.073 0.120
	average st. dev.		0.068 0.016	0.130 0.053

ACKNOWLEDGEMENTS

This work was not possible without the help of V. Jeammet, curator at the Greek, Etruscan and Roman Antiquities in Louvre Museum, general curator of the exhibition, and her co-workers J. Becq and N. Mathieux, who kindly provided figurines to be sampled.

We are indebted to G. Querré, the former leader of the thermoluminescence group, who initiated the luminescence dating of the Tanagras and introduced optical dating to the laboratory.

We would like to thank our colleagues from the laboratory involved in the Tanagra project: A. Bouquillon (clay studies), S. Collinart, L. Troalen, E. Van Elsande (pigment analyses), D. Bagault (photography), and M. Dubus (Xrays diffractometry) for discussion and comments. We are grateful to A. Bouquillon, leader of the ceramic-dating group, for support of the luminescence unit.


U-Chan Chung Seu

Catherine Elissalde

Stéphane Mornet

Mario Maglione

Controlling internal barrier in low loss BaTiO3 supercapacitors

Frittage Flash, PNF2-CNRS, 118 route de Narbonne, 31062 Toulouse cedex 09, France Supercapacitor behavior has been reported in a number of oxides including reduced BaTiO 3 ferroelectric ceramics. These so-called giant properties are however not easily controlled. We show here that the continuous coating of individual BaTiO 3 grains by a silica shell in combination with Spark Plasma Sintering is a way to process bulk composites having supercapacitor features with low dielectric losses and temperature stability. The silica shell acts both as an oxidation barrier during the processing and as a dielectric barrier in the final composite.

Commercial supercapacitors are using polycrystalline materials in which metallic grains or foils (like tantalum) with extended specific surface are covered by their respective oxides (like tantalum oxide). The resulting barrier layers locate electronic and ionic space charges, increasing the capacitance density. While such barrier layer supercapacitors are of everyday use, 1,2 an improved control of the interfaces barrier in such advanced materials would lead to increased functionalities: higher charge storage ability, lower losses and increased cut-off frequency. We show here that individual coating of ferroelectric BaTiO 3 grains by a silica shell followed by a fast sintering in reducing conditions leads to such stable supercapacitor behavior in electroceramics. The tuning of a dielectric shell at the grain scale, as compared with standard routes, offers several opportunities for electroceramics.

Actually, there is an intense wave of research on materials displaying extreme dielectric parameters. Many compounds such as CaCu 3 Ti 4 O 12 and LuFe 2 O 4 have been quoted

as "giant" permittivity materials. Well known and widely used ceramics like BaTiO 3 show huge effective dielectric parameters when appropriately substituted or when decreasing the grain size under fast sintering conditions. [3][4][5][6] In all cases, the apparent dielectric permittivity is high at room temperature (> 10 5 in some cases) and is temperature independent in a broad range before falling down at low temperatures. While these features are very appealing for applications, some parameters are still not controlled: the dielectric losses are high and low frequency contributions to the dielectric permittivity are ascribed to dc conductivity. In CaCu 3 Ti 4 O 12 ceramics, this dielectric behavior was ascribed to an internal barrier layer (IBLC) at the grain boundaries. 7 Within this model, the measured huge capacitances are not increased by intrinsic dielectric parameters but by geometrical confinement of the dielectric at the grain boundaries. All the above mentioned materials can be described using this IBLC model, each of them having different free charges able to increase the inner grain conductivity. [START_REF] Maglione | Charge Transfer and Vibronic States in Ionic-Covalent Systems" Theory, Experiment, and Applications[END_REF] However, the key parameter which is the grain boundary interface, is still not well understood and controlled. This is a strong drawback in order to reproduce the dielectric performances and to reduce the dielectric losses of IBLC ceramics. We propose here to solve this issue by using advanced processing of core@shell particles sintered using Spark Plasma Sintering (SPS). 9,10 By means of homogeneously coated BaTiO 3 grains with an amorphous silica shell, we were able to obtain composite ceramics of apparent permittivity of 2.10 5 and losses of only a few percent at 10 4 Hz which are stable against temperature.

The two main conditions for the internal barrier layer to increase the overall dielectric permittivity of oxides are the inner grain conductivity and the insulating grain boundary. The former is related to the amount of charged defects which are intentionally introduced prior to the precursor reaction or during the sintering step under reducing conditions. 6,11 The stabilization of such charged lattice defects in a reproducible way in BaTiO 3 ceramics, is still a challenge. Even the way to compensate for heterovalent substitution like La 3+ on the Ba 2+ site of BaTiO 3 is still a matter of debate. 11,12 The control of the grain boundaries (i.e. tuning their morphologies and electrical properties) is also not an easy task. Indeed, depending on the size and the chemical state of the starting grains, the sintering step of the bulk ceramics can strongly alter shape, charges and electronic state of the grain boundaries. For BaTiO 3 , many discrepancies have been quoted in the literature depending on the processing route. [13][14][15] This is even more true when fast sintering processes under reducing atmosphere like Spark Plasma Sintering are used. This is why we decided to coat each BaTiO 3 grain prior to the sintering step in order to improve the reproducibility of the grain boundaries in dense ceramics. The first step is to produce core@shell powders made of BaTiO 3 cores surrounded by a continuous and homogeneous shell of amorphous silica (SiO 2 ). We previously showed that the thickness of the shell can be tuned continuously from few nanometers up to tenths of microns. This core@shell architecture was efficient to decrease the dielectric losses of the final ceramics and also to achieve multifunctionality on including nanoparticles in the coating. 9, 16 For the present purpose of achieving extremely high and reproducible dielectric permittivity ceramics, we started from standard BaTiO 3 particles of 500 nm diameter purchased from Sakai Chemical Co. (Japan). All the particles have been individually coated with a homogeneous amorphous silica shell of 5 nm thickness using a method derived from the Stöber process, as described in detail elsewhere. 9 Core@shell BaTiO 3 (500 nm)@SiO 2 (5 nm) particles are denoted as

[BT500@Si5].
Once the starting core@shell particles are available, the sintering step should also be optimized since any porosity is detrimental for getting large capacitance density. Improved two-step sintering has been used but long dwelling times are required (950°C, 12h) and the resulting density does not reach 85%. 17 As to build an IBLC ceramic, we then need to induce charged point defects in the ceramics. It has been previously shown that SPS, which is operated under reducing atmosphere, is able to create such defects as Ti 3+ associated to oxygen vacancies. 6 The SPS experiments were performed in a Dr Sinter SPS-2080 Syntex INC Japan. A heating rate of 100°C.min -1 was used to reach the final temperature of 1100°C under vacuum. A uniaxial pressure of 50 MPa was applied immediately before and until completion of the temperature rising step. The as obtained ceramics have a surface layer contaminated by carbon coming from the graphite die in which the pellets are pressed. A postannealing step at 800°C for 12h in air is performed for decarbonation of the samples; the density reached 97% for all of them. After this first post-annealing treatment, the pellet remains navy blue stemming from the reduction state of titanium (this color being usually ascribed to Ti 3+ -VO centers).

A second identical thermal treatment did not allow re-oxidizing the sample, highlighting the efficiency of silica as a barrier against re-oxidation. Considering the extreme SPS sintering conditions, the silica in the final pellet, can be viewed as a glassy matrix in which the BT grains are embedded as evidenced by High Resolution SEM. The grain size remains close to the initial one (i.e. 500 nm) and the fused aspect on the grains surface is attributed to the glassy silica matrix (Fig. 1).

Having confirmed that charged defects are indeed stable in the ceramics, we next need to check the silica shell at the grain boundaries. In a separate paper, we used in situ Transmission Electron Microscopy (TEM), High Resolution TEM and XPS, showing the preservation of the silica coating under the mentioned SPS conditions [submitted]. We point out that cross diffusion between silica and BaTiO 3 is moderate and that amorphous silica is still surrounding the grains in the ceramics. The capacitance and dielectric losses of the samples were recorded as a function of temperature (4 K-450 K) and frequency (100 Hz-10 MHz). On figure 2, we report the real part of the dielectric permittivity (') and the dielectric losses (tg ) in the temperature range 300-450 K for several spot frequencies. These parameters are only apparent ones and even if they should be considered with some care, because of the IBLC model, we will, however, discuss them for sake of comparison with the existing literature. First, the room temperature dielectric permittivity is in the range of the so called giant dielectric materials (~ 2.10 5 at f= 10 4 Hz) meaning that our core shell particles are indeed leading to IBLC ceramics (Fig. 2a). The main achievement here is that the dielectric losses are much smaller than the ones previously reported i.e. tg ~ 5% at f= 10 4 Hz instead of 100% in the existing literature (Fig. 2b). Moreover, the absence of apparent conductivity contribution at low frequencies shows a small conductivity across the whole ceramic even at T>100ºC. These two main features, lower dielectric losses and the absence of a conductivity contribution may be ascribed to the silica coating. Another interesting feature is the small maximum in the dielectric permittivity at the ferroelectric/paraelectric transition temperature T= 400 K (Fig. 2a). This means that the bulk ferroelectric properties of the material are not totally cancelled by the IBLC effect.

To get a more quantitative insight into this IBLC contribution to the increased dielectric parameters, one should follow the dielectric relaxation which inevitably occurs for all of similar materials. [START_REF] Maglione | Charge Transfer and Vibronic States in Ionic-Covalent Systems" Theory, Experiment, and Applications[END_REF] Such relaxation clearly appears on figure 3a where the real part of the dielectric permittivity bends down sharply at low temperatures. As for every supercapacitor, the high frequency permittivity is thus much lower than the low frequency one. The relaxation is better seen on the losses which undergo a strong maximum and whose temperature is shifted with the operating frequency (Fig. 3b). Both the shape and the shift of this maximum stem from a dielectric relaxation which is nothing but the relaxation of the space charges located at the grain boundaries. The activation energy of this space charge may be mapped out on plotting the characteristic frequencies and temperatures in an Arrhenius plot (Fig. 3a inset).

The linear evolution evidences the exponential behavior of the relaxation frequency

        kT E f f exp 0
where f 0 is the high temperature extrapolation of the relaxation frequency and E its activation energy; E= 0.15 eV which agrees with similar relaxations observed in many other materials. [START_REF] Maglione | Charge Transfer and Vibronic States in Ionic-Covalent Systems" Theory, Experiment, and Applications[END_REF] This supports space charge relaxation at the IBLC as a probable mechanism. Indeed, 0.15 eV is in the range of the Ti 3+ level position below the conduction band of BaTiO 3 . Electron hopping between such defects and their localization at the interface is thus a way to artificially increase the dielectric properties. 18,19 Compared with IBLC ceramics referred in literature, our core@shell composites exhibit unprecedented behavior as they combine both giant permittivity and stable dielectric losses as low as 5% at 10 kHz. In a similar way, our structures can be called ISBLC (S for silica). The occurrence of the Curie temperature maximum indicates that the bulk properties of the composite are not cancelled by the silica barrier layer. Microstructural investigations allowed to describe the final composite structure as ferroelectric grains embedded in a glassy matrix. SPS does not only help for densification but allows to quench the amorphous silica in a glassy state thanks to the extremely fast sintering and cooling rates. In turn, this glassy silica cage allows an irreversible reduction of the BT core (induced by the reducing SPS sintering conditions). As a result, the SPS allows acting on both the core and the shell. The silica is both a dielectric barrier and a re-oxidation barrier leading to low loss supercapacitor. In addition, its thickness can be easily controlled leading to reproducibility and reliability of our composite properties. The supercapacitor behavior has been ascribed to space charge relaxation associated to Ti 4+ reduction.

1

  Controlling internal barrier in low lossBaTiO 3 supercapacitors U-C. Chung, C. Elissalde, S. Mornet and M. Maglione ICMCB -CNRS, Université Bordeaux, 87 Avenue du Dr A. Schweitzer, F-33608 Pessac cedex, France C. Estournès CNRS; Institut Carnot Cirimat; F-31062 Toulouse, France et Plateforme Nationale de

FIG. 1 .

 1 FIGURE CAPTIONS

Fig. 1

 1 Fig. 1

Fig. 2

 2 Fig. 2

Fig. 3

 3 Fig.3