Oscar Acosta

Jurgen Fripp

Vincent Doré

Pierrick Bourgeat

Jean-Marie Favreau

Gaël Chételat

Andrea Rueda

Victor L Villemagne

Vincent Dore

Cassandra Szoeke

David Ames

Kathryn A Ellis

Ralph N Martins

Colin L Masters

Christopher C Rowe

Erik Bonner

Florence Gris

Di Xiao

Parnesh Raniga

Vincent Barra

Olivier Salvado

Cortical surface mapping using

Keywords: Cortical mapping, Surface registration, Partially inflated surfaces, Cortical thickness, Brain Atrophy, Shape context, Statistical analysis of populations, Alzheimer's disease, Freesurfer

Introduction

Non--invasive imaging techniques such as magnetic resonance imaging (MRI) allow the investigation of many of the functional and morphological changes in the brain occurring with age or during pathological processes. However, reliable statistical analysis of available data over a large number of individuals are required before detecting subtle differences over time or between individuals and groups. The distribution differences of specific tissues such as white matter (WM), gray matter (GM) or cerebrospinal fluid (CSF), volume changes, or any other signal coming from imaging devices or intermediate pre--processing stages have been used in the past, unravelling some of the complex pathological changes. Statistical analysis per region (Lerch et al. 2005), voxels [START_REF] Ashburner | Voxel--based morphometry-the methods[END_REF], (Ziolko and textitet al 2006) or vertices in a surface representation of the brain (A.M. Dale, B. Fischl, and M.I. Sereno 1999), (B. Fischl, M. I Sereno, and A. M Dale 1999) can be used to assess the differences between groups using, for instance, a general linear model (GLM) to the available data (intensity distribution, deformation, etc.) from all subjects at each region (voxel or vertex) with different covariates such as age, years of education, gender, diagnosis, cognitive scores, etc. Some of these techniques have quantitatively demonstrated different biomarkers or predictors for Alzheimer's Disease (AD), showing a specific pattern of atrophy [START_REF] Hua | 3D characterization of brain atrophy in Alzheimer's disease and mild cognitive impairment using tensor--based morphometry[END_REF], [START_REF] Querbes | Early diagnosis of Alzheimer's disease using cortical thickness: impact of cognitive reserve[END_REF], hypometabolism [START_REF] Chételat | Direct voxel--based comparison between grey matter hypometabolism and atrophy in Alzheimer's disease[END_REF]), β --amyloid load [START_REF] Mikhno | Voxel--based analysis of 11C--PIB scans for diagnosing Alzheimer's disease[END_REF] or combined hypometabolism and atrophy [START_REF] Kawachi | Comparison of the diagnostic performance of FDG--PET and VBM--MRI in very mild Alzheimer's disease[END_REF], [START_REF] Villain | Relationships between Hippocampal Atrophy, White Matter Disruption, and Gray Matter Hypometabolism in Alzheimer's Disease[END_REF]) and sulcal modifications [START_REF] Zhan | An entorhinal cortex sulcal pattern is associated with Alzheimer's disease[END_REF]) using MRI or Positron Emission Tomography (PET).

Different software tools, grouping a wide range of methods to perform studies of cortical and subcortical anatomy are nowadays available. Most of them are also able to carry out further population analysis. Statistical Parametric Mapping (SPM) 1 or FreeSurfer (FS) 2 represent two of those automated tools largely referenced in the literature. Although they share the same goals, the comparisons are performed in different ways. SPM is a voxel--based approach in which the images are realigned, spatially normalised into a standard space, and smoothed before performing statistical analysis at a voxel level [START_REF] Ashburner | Voxel--based morphometry-the methods[END_REF]. In contrast, FreeSurfer is a surface--based approach performing first the reconstruction of the brain cortical surface from structural MRI data (B. Fischl, M. I Sereno, and A. M Dale 1999).

FreeSurfer provides many anatomical analysis tools, including the representation of the cortical surface, segmentation of the brain, skull stripping, bias field correction, nonlinear registration of the cortical surfaces into a standard space, labelling of regions and statistical analysis of populations. Comparison of those different tools has been performed in the past for specific tasks such as segmentation [START_REF] Klauschen | Evaluation of Automated Brain MR Image Segmentation and Volumetry Methods[END_REF]. Overall, they face similar methodological challenges for interindividual comparisons.

The primary issue for interindividual comparison is the mapping or registration of a given population to a common space. This is a particularly challenging task because of the brain convoluted geometry and the high interindividual variability (J.-F. [START_REF] Mangin | A framework to study the cortical folding patterns[END_REF]). Various registration approaches have been proposed, including intensity--based in the voxel space [START_REF] Rueckert | Nonrigid Registration Using Free--Form Deformations: Application to Breast MR Images[END_REF]), (T. [START_REF] Vercauteren | Non--parametric Diffeomorphic image Resgistration with the Demons Algorithm[END_REF] or shape--based in a surface representation of the brain [START_REF] Yeo | Spherical Demons: Fast Surface Registration[END_REF][START_REF] Yeo | Spherical Demons: Fast Diffeomorphic Landmark--Free Surface Registration Medical Imaging[END_REF], [START_REF] Eckstein | Generalized Surface Flows for Deformable Registration and Cortical Matching[END_REF]). Although intensity--based non--rigid registration methods work well in intrasubject registration, cortical folds do not match well and may not effectively address the issues arising from the variability across a population [START_REF] Hellier | Retrospective evaluation of intersubject brain registration[END_REF]. Alternatively, cortical surfaces registration allows anatomically meaningful features to constrain the transformations (Tosun et al. 2004), which improves the reliability of the registration and rather favors the alignment of functional regions. Although variability in high level foldings, gyri and sulci, exists, the cortical shapes tend to present similar patterns at coarse levels [START_REF] Tosun | Mapping techniques for aligning sulci across multiple brains[END_REF]. A good matching may be consequently found by the alignment of the lobes and major folding patterns at a coarse level across the brains, thereby registering cortical surfaces instead of volumetric data.

Reliable cortical surface registration for statistical comparisons raises different methodological questions which constrain any pre--processing step. The main issues are related with i) the preservation of gyri and sulci and with ii) the topology correction. Firstly, since the folding patterns correlate well with functional and anatomical regions between individuals [START_REF] Fischl | Cortical folding patterns and predicting cytoarchictecture[END_REF]) a reliable detection of gyri and sulci during the segmentation step is required. Various methods have been proposed: Hutton (Hutton et al. 2008) used a layering method in the voxel space based on mathematical morphology to detect deep sulci. In [START_REF] Acosta | Automatic delineation of sulci and improved partial volume classification for accurate 3D voxel--based cortical thickness estimation from MR[END_REF], a method is proposed to improve sulci detection by using a distance based cost function from WM in a post processing step, but no verification of topology is done.

In [START_REF] Rueda | Topology--corrected segmentation and local intensity estimates for improved partial volume classification of brain cortex in MRI[END_REF] the approach is improved by introducing topology constraints but also in a post--processing step, obtaining both reliable pure tissue and partial volume images segmentations. In [START_REF] Cardoso | LoAd: A locally adaptive cortical segmentation algorithm[END_REF] an explicit model of partial volume classes is introduced within the segmentation step and a locally varying MRF--based model is used to locally modify the priors for enhancement of gyri, but the topology is not guaranteed yet. Secondly, the GM can be considered as a folded sheet built upon the WM and it is often assumed that if the midline hemispheric connections were artificially removed, the cortex would have the topology of a hollow sphere (neither handles nor tunnels). This assumption must be preserved throughout all the processing pipeline. Topological operators and constraints have to be used to correct and achieve accurate cortical tissue representations either in the voxel space during the segmentation or in the surface representations [START_REF] Ségonne | Geometrically Accurate Topology--Correction of Cortical Surfaces Using Nonseparating Loops[END_REF], [START_REF] Ségonne | Active Contours Under Topology Control --Genus Preserving Level Sets[END_REF], [START_REF] Bazin | Topology Correction Using Fast Marching Methods and Its Application to Brain Segmentation[END_REF], [START_REF] Han | Topology correction in brain cortex segmentation using a multiscale, graph--based algorithm[END_REF], [START_REF] Kriegeskorte | An Efficient Algorithm for Topologically Correct Segmentation of the Cortical Sheet in Anatomical MR Volumes[END_REF]. The main drawback of correcting topology in the voxel space lies in the modifications introduced by adding or removing voxels compared to the size of the structures, besides the non unicity of the handle definition. Some other approaches used a combined voxel--based strategy [START_REF] Jaume | Open Topology: A Toolkit for Brain Isosurface Correction[END_REF], [START_REF] Zhou | Topology Repair of Solid Models Using Skeletons[END_REF] and, as we mentioned before, other operates directly with the surfaces (CLASP (J. S [START_REF] Kim | Automated 3D extraction and evaluation of the inner and outer cortical surfaces using a Laplacian map and partial volume effect classification[END_REF], BrainVISA (J-F. [START_REF] Mangin | From 3D magnetic resonance images to structural representations of the cortex topography using topology preserving deformations[END_REF], so as Freesurfer (A.M. Dale, B. Fischl, and M.I. Sereno 1999), (B. Fischl, M. I Sereno, and A. M Dale 1999), (B. Fischl and A. M Dale 2000)). FreeSurfer incorporates mechanisms to prevent self-intersection of surfaces or topology correction, additionally imposing smoothness constraints, but at expense of computational cost.

The goal of topology correction in surface--based methods is the elimination of defects appearing in the mesh such as tunnels and handles after the mesh generation (B. Fischl, M. I Sereno, and A. M Dale 1999), (Florent Ségonne and Bruce Fischl 2007), [START_REF] Jaume | Open Topology: A Toolkit for Brain Isosurface Correction[END_REF], [START_REF] Zhou | Topology Repair of Solid Models Using Skeletons[END_REF], [START_REF] Hong | Conformal virtual colon flattening[END_REF]. Some of them inflate the surface to detect the topological incoherences (B. Fischl, M. I Sereno, and A. M Dale 1999), (Florent Ségonne and Bruce Fischl 2007). Nevertheless, the optimization step for spherical inflating is non--deterministic and its complexity cannot be evaluated. Moreover, the corrections are not optimal from the topological point of view as the surfaces are not corrected according to the shortest loops that corresponds to the handles and tunnels. In this paper we addressed the issue of topology correction with a method based on the detection and optimal cutting of non-separating loops in the mesh.

Concerning the surface based non--rigid registration, some of the previously proposed approaches have used a parametric representation on a sphere, obtained either through iterative relaxation like in Freesurfer (B. Fischl, M. I Sereno, and A. M Dale 1999) or conformal mapping [START_REF] Angenent | On the Laplace--Beltrami operator and brain surface flattening[END_REF], [START_REF] Gu | Genus zero surface conformal mapping and its application to brain surface mapping[END_REF]). In the spherical domain, the registration strategy relies on features such as sulcal landmarks or mean curvature, or convexity [START_REF] Yeo | Spherical Demons: Fast Diffeomorphic Landmark--Free Surface Registration Medical Imaging[END_REF], [START_REF] Yeo | Spherical Demons: Fast Surface Registration[END_REF] , [START_REF] Ziolko | Evaluation of voxel--based methods for the statistical analysis of PIB PET amyloid imaging studies in Alzheimer's disease[END_REF]. In those approaches the global shape information is not explicit and the registration is dependent on the selected features. [START_REF] Shi | Direct cortical mapping via solving partial differential equations on implicit surfaces[END_REF] proposed the computation of a direct map from the source to the target constrained by the sulcal landmark curves avoiding any parameterization. It assumes nevertheless that the sulcal landmarks curves are known and can be matched, which is not always the case. An intermediate unfolded (inflated) representation of the brain provides an alternative to expose hidden sulci and to simplify the geometry for cortical mapping preserving its global shape at a coarse level such as the main lobes and folds. Known as partially flattened surfaces (PFS) [START_REF] Drury | Computerized Mappings of the Cerebral Cortex: A Multiresolution Flattening Method and a Surface--Based Coordinate System[END_REF] , (B. Fischl, M. I Sereno, and A. M Dale 1999), [START_REF] Pons | Area Preserving Cortex Unfolding[END_REF] they can also provide an intermediate step to further conformal mapping [START_REF] Tosun | Mapping techniques for aligning sulci across multiple brains[END_REF]. In a previous work we compared different methods for obtaining PFS in terms of their ability to preserve areas and angles [START_REF] Bonner | A quantitative comparison of three methods for inflating cortical meshes[END_REF]. Compared to the original brain surface representations, the PFS simplifies the determination of shape correspondences between subjects while avoiding the spherical mapping procedure [START_REF] Eckstein | Generalized Surface Flows for Deformable Registration and Cortical Matching[END_REF]. A meaningful descriptor of shape in terms of brain anatomy may be used, thereby, as a similarity criterion.

In this paper we present a new cortical surface processing pipeline for vertex--wise inter-individual analysis. We provide a comprehensive approach which performs the surface cortical mapping steps through MRI segmentation, surface generation, topology correction, inflation and surface matching, yielding a representation of a whole population in a common space. The matching between individuals' flattened surfaces (PFS) is based on the highly discriminative scale--invariant properties of the shape context (SC) (S. [START_REF] Belongie | Matching with shape contexts. Content--based Access of Image and Video Libraries[END_REF], previously explored in [START_REF] Acosta | 3D shape context surface registration for cortical mapping Proc[END_REF]. The similarity between points is a combination of shape context and the information of sulci automatically extracted during an inflation step (sulci depth). The interpolation is obtained with the Thin Plate Spline (TPS) [START_REF] Bookstein | Principal warps: thin--plate splines and the decomposition of deformations[END_REF]). The major contribution of the shape context lies in the global shape characterization at a local level for each single point. Our method produces an anatomically meaningful and scale--invariant matching between the lobes and the major folding patterns.

The C++ code and binaries for segmentation, cortical thickness estimation (CTE)

proposed in [START_REF] Acosta | Automated voxel--based 3D cortical thickness measurement in a combined Lagrangian--Eulerian PDE approach using partial volume maps[END_REF] and the topology correction are fully available at our CSIRO software website3 . The proposed pipeline, allows also the comparison of different imaging biomarkers for other neurodegenerative diseases. In the remainder of the paper we explain the following steps: generation of genus zero surfaces (Section 2.1), PFS (Section 2.2) and a description of the shape context (Section 2.3.1). The validation experiments using real MRI are then described and a clinical application in an Alzheimer's disease study is also presented (Section 3.4) where we tested the ability of cortical thickness as a biomarker for atrophy to compare two populations.

[ Figure 1 : Overall process for cortical mapping]

Materials and Methods

The proposed method consists of several stages as summarized in Fig. 1: Firstly, 3D T1-weighted MR images are classified into GM, WM and CSF in their original space using an expectation maximisation segmentation (EMS) algorithm [START_REF] Van Leemput | Automated model--based bias field correction of MR images of the brain[END_REF], as described in [START_REF] Bourgeat | AB burden in the temporal neocortex is related to hippocampal atrophy in non--demented elderly[END_REF]). The EMS computes probability maps for each tissue type, which are then discretized by assigning each voxel to its most likely tissue type. This scheme has been improved in [START_REF] Bourgeat | AB burden in the temporal neocortex is related to hippocampal atrophy in non--demented elderly[END_REF]) allowing the use of multiple (n) atlases for subsequent segmentations and parcellations, in order to improve the reliability of the priors. Thus, each individual scan is segmented n times, once for each of the n atlases. The resulting n segmentations per patient are then combined using a voting scheme to provide consensus segmentation according to [START_REF] Aljabar | Multi--atlas based segmentation of brain images: Atlas selection and its effect on accuracy[END_REF]. Similarly, for each patient, the n propagated Automated Anatomical Labeling (AAL) [START_REF] Tzourio--Mazoyer | Automated anatomical labeling of activations in SPM using a macroscopic anatomical parcellation of the MNI MRI single--subject brain[END_REF] and the Internet Brain Segmentation Repository (IBSR) labeling 4 , are also combined using the same voting scheme.

This step ensures that outliers resulting from poor spatial normalization are excluded, and only regions of high confidence are preserved. Secondly, based on these pure tissue segmentations, a further maximum a posteriori classification of voxels into pure tissues WM, GM and CSF and mixed tissues WM/GM and GM/CSF along the previously computed GM interface is performed, which results in a GM partial--volume coefficients (GMPVC) image used for further cortical thickness estimation [START_REF] Acosta | Automated voxel--based 3D cortical thickness measurement in a combined Lagrangian--Eulerian PDE approach using partial volume maps[END_REF]). Topology--constraints are introduced in the voxel classification assuming that the GM is a continuous layer covering a WM, homotopic to a filled sphere as presented previously in [START_REF] Rueda | Partial Volume Estimation of Brain Cortex from MRI Using Topology--Corrected Segmentation[END_REF]. A topology preserving dilation of the WM over the GM improves the robustness of delineation of mixed GM/CSF voxels in deep sulci. The propagated IBSR labeling is then used as a mask to separate left and right hemispheres, after which the generation of 3D polygon meshes representing the WM/GM interface is performed with the marching cubes algorithm [START_REF] Lorensen | Marching cubes: A high resolution 3D surface construction algorithm[END_REF]. Tunnels and handles resulting from the meshing are removed with a method based on the detection and correction of non-separating loops (Fig. 3). In order to reduce the surface convolution and the inter--individual variability, the surfaces are then inflated to a similar shape (Fig. 4). The inflated surfaces, referred to as PFS, are finally elastically registered towards the common template using the shape context--based non rigid registration (Fig. 5). The registration imposes a standardized coordinate system on both surfaces, allowing the vertex--wise comparison of the whole population. In this paper, we detail and validate the steps 5 to 7 in Fig. 1 and illustrate the results with an example of clinical application in which the cortical thickness was computed using our voxel--based approach (Acosta et al. 2009) propagated to the registered meshes. We compared the differences in cortical thickness between a population of Alzheimer's disease (AD) patients and a control group of healthy elderly individuals (HC).

Correction of topology and generation of genus zero surfaces

After the surfaces are generated by the Marching Cubes algorithm [START_REF] Lorensen | Marching cubes: A high resolution 3D surface construction algorithm[END_REF] the topology is corrected to meet the assumed topological properties of the brain, homeomorphic to a sphere. Handles and tunnels appear as topological defects to be removed.

To this end, we implemented an iterative topological correction inspired by (Erickson and Hard--Peled 2002) allowing to remove one by one the unwanted handles and tunnels.

Let g be the genus of the surface M , defined as the number of handles and tunnels in M . g can be easily computed using the Euler's characteristic

) M ( χ as 2 ) M ( 2 χ - = g
, where

) M ( ) M ( ) M ( ) M ( T E V + - = χ and ) M ( V , ) M ( E and ) M ( T
are respectively the number of vertices, edges and triangles in M .

The algorithm (Algorithm 1) iteratively corrects the topology, reducing the genus of M by one at each step. After g corrections, filling the tunnels and cutting the handles, the new surface is homeomorphic to a sphere.

Observe that in a first step the global algorithm is applied aimed at closing all the possible boundaries of the original surface. The removal of a handle or a tunnel is performed by first computing an approximation of the shortest non--separating loop at each topological defect (detailed below and summarized in the Algorithm 2). A non--separating loop is a connected path that does not divide the surface into two connected components. The surface M is cut along the n edges of this identified loop, therefore adding n edges and n vertices to the mesh. At this step, two new boundaries are created but the Euler's characteristic is not modified. The final correction is performed by sealing the created boundaries with small discs, thus adding n edges n triangles and one vertex (we add an umbrella). This increases the Euler's characteristic by one, globally reducing the genus by one.

Algorithm 2 describes the computation of the shortest non--separating loop. We first compute a set of basepoints B by propagating a wavefront from a randomly selected vertex (a seed point) in M . B is defined as the set of vertices where the two wavefront boundaries meet (Figure 2a). The complement of B in M is a disc, implying that any non--separating loop on M contains at least one point of B . Figure 2b is an example of cut locus on a 2--genus surface.

Using this property we compute the non--separating loops from all the vertices contained in the path B . This has been proved by (Erickson and Hard--Peled 2002) as a good approximation of the shortest non--separating loop in M . Unlike the original approach, in which each point of B is used to compute a non--separating loop, we perform this computation using propagations from a non--trivial set of points, which speeds--up the computation. As illustrated in Figure 2b, the cut locus B is subdivided into a finite (and small) number of paths, like B S in Figure 2c. For each of these paths, we compute a wavefront from the contained points into the complement of B in M . When the two boundaries of this wavefront meet, an associated loop is computed from the junction point going back to the initial wavefront seed.

This, at each side of the junction. The non--separating property of this loop is verified by checking the number of connected components of the complement. Selecting a good approximation of the shortest non--separating loop at each step guarantees a minimal structural modification of the surface M . Figure 3 shows topological corrections applied to a tunnel (figure 3a, 3b) and to two handles (figure 3c,3d).

[Figure 2]

[Figure 3]

Partially flattened surfaces

A number of methods have been proposed in the literature for unfolding or flattening cortical surfaces [START_REF] Drury | Computerized Mappings of the Cerebral Cortex: A Multiresolution Flattening Method and a Surface--Based Coordinate System[END_REF], (B. Fischl, M. I Sereno, and A. M Dale 1999) , [START_REF] Pons | Area Preserving Cortex Unfolding[END_REF]. In order for a flattening to be useful, it must preserve local and global metrics such as triangle angles and areas. For our purposes we implemented a method based upon CARET [START_REF] Drury | Computerized Mappings of the Cerebral Cortex: A Multiresolution Flattening Method and a Surface--Based Coordinate System[END_REF], [START_REF] Tosun | Mapping techniques for aligning sulci across multiple brains[END_REF]. This method was compared with the method implemented in Freesurfer [START_REF] Bonner | A quantitative comparison of three methods for inflating cortical meshes[END_REF]. In our implementation, a cortical surface is iteratively deformed at each vertex according to 1, where t i m is the position of vertex i for iteration t , λ is a scalar in the range [0, 1] and t i m represents the average vertex position of t i m , given by

t i t i t i m m m λ λ + - = + ) 1 (
t j i N j t j t j t i c S S = m ∑ ∑ ∈ 1 (1)
where i N is the set of all triangles containing i m , j c is the center of triangle j and j S is its area. Typically a factor of 0.9 is chosen for λ , giving more weigh to the averaged vertex position. This process moves all mesh vertices towards the weighted average of the centers of their surrounding triangles. The deformation progresses until the global mean curvature k drops below a predefined threshold [START_REF] Tosun | Mapping techniques for aligning sulci across multiple brains[END_REF] or in our case using a predefined number of iterations (we used typically 400), obtaining similar shapes between moving and fixed surfaces.

∑ = + - = ] ,..., 1 [ 1 ] [ N t t i t i i m m depth ( 2
)
where N is the number of iterations. [Figure 4]

Surface--based non--rigid registration

Surface based non--rigid registration methods face several challenges including: i) the choice of similarity criterion and ii) the matching and global optimization procedure [START_REF] Audette | An algorithmic overview of surface registration techniques for medical imaging[END_REF]. The first one refers to the type of information extracted from the 3D surface, namely the description of local or global shape to represent the similarity. The latter challenge concerns the exploitation of the similarity information to find the best matching between the two surfaces. In this context, the goal of the registration is to determine the transformation such that for a finite set of control points, any control point of a moving surface M , is mapped onto the corresponding control point of a fixed surface M .

[Figure 5] Fig. 5 illustrates the method used to register the moving

p i i ... 1 } m { M = =
and fixed

q j j ... 1 } f { F = =
PFS surfaces, respectively. F may be a template towards which we register the whole population. Let

u k k ... 1 s } m { M = = and v l l ... 1 } f { F = =
be the corresponding simplified surfaces, represented in two subsets of characteristic points (control points) of M and F , such that M M s ⊆ , F F s ⊆ . After establishing bijective correspondences between the sets s M and s F , exploiting a given similarity metric (in this paper the shape context), the mapping is computed as a set of transformations

u k k ... 1 } T { T = = ) R R (T 3 3 k → ∈ such that k k m , m T f ∀ =  k k .

Shape context

The shape context is a shape descriptor that, for a single surface, captures the distribution of points over relative positions of the global shape points. This characterization is invariant to scale and rigid transformations, naturally leading to a highly discriminative and robust score for measuring shape similarities. The shape context was first introduced in (S. [START_REF] Belongie | Matching with shape contexts. Content--based Access of Image and Video Libraries[END_REF] within the 2D pattern recognition field, and aimed to match point clouds representing similar patterns. Further modifications to the method appeared in (G.

Mori, S. Belongie, and J. Malik 2005) and 3D matching of features have been extended to work with thoracic images [START_REF] Urschler | Registering 3D Lung Surfaces Using the Shape Context Approach[END_REF]. A first implementation of 3D Shape Context for brain registration was presented in [START_REF] Acosta | 3D shape context surface registration for cortical mapping Proc[END_REF].

Let s M and s F the two simplified PFS surfaces corresponding to the meshes M and F to be registered. s M and s F constitute the u control points (typically u = 3000 points) of the moving and fixed surfaces, respectively, to be matched and for which the shape context is computed (Fig. 5). Here, they are obtained using an algorithm that iteratively contracts vertex pairs while minimizing geometric errors [START_REF] Garland | Surface simplification using quadric error metrics[END_REF]. For a given point

s ... 1 M } m { ∈ = u k k , its shape context is the 3D histogram of the relative 3D polar coordinates ) , , ( φ θ r 
of the remaining 1 u points. As in [START_REF] Belongie | Shape matching and object recognition using shape contexts[END_REF], in order to be more sensitive to nearby points, we use a log--polar coordinate system. In our case we build a 3D histogram with r R  equally spaced log--radius bins and θ u and ϕ v equally spaced angle bins.

Since the intracranial volume may considerably differ across a population, an additional normalization of the shape context is performed to obtain scale invariance. To this end, all the

radial distances } { , ) m , m ( ... 1 , l k l k dist u l k ≠ =
, are normalized by the mean distance between all the point pairs in the shape:

∑ ≠l k l k dist u , l k ) m , m ( 1
(3)

Cost function and matching

The shape context cost SC j i, C of matching the point i m from the simplified surface s M

with the point j f from s F is given by the 2

χ statistic as [ ] ∑ - K = k j i j i SC j i, (k) g + (k) h (k) g (k) h = C 1 2 2 1 (4) where ) (k h i and ) (k g j
are the histograms (shape context) at the points i and j of the moving s M and fixed s F shapes, respectively, and K is the total number of bins.

In addition to the description given by the shape context, we use the sulcal depth map computed at each point according to equation ( 2) to enrich the matching between the points with local information about sulci and gyri. Thus, in a second iteration we recompute the depth map cost function SD j i, C as :

| | N r s = C j j SD j i, - (5)
Where i s and j r are the sulcal depths at points i and j of the moving s M and fixed s F shapes respectively. Both are normalized by the brain size. N is the highest difference between the floating and the fixed normalized sulcal depths; this enables the comparison between the sulcal depth and the shape context cost. The global cost function becomes then:

SD j i, SC j i, j i, βC + C α = C (6)
where . This is done within a one to one point matching step with the Hungarian algorithm [START_REF] Kuhn | The Hungarian method for the assignment problem[END_REF]. After the correspondences are found, the set of

transformations u k ,.., 1 k } {T T = =
for each point are computed. Finally, a transformation using a thin plate spline model (TPS) [START_REF] Bookstein | Principal warps: thin--plate splines and the decomposition of deformations[END_REF]) calculated from the control points, achieves the interpolation of the corresponding moving shape M onto F . In our implementation, the number of points in s M is lower than in s F to avoid miscorrespondences or topological errors due to crossing points between s M and s F . Fig. 6 shows an example of labeling, propagated from a template to a single individual's surface after registration of the corresponding PFS.

Experiments and results

This section describes the experiments performed to validate each step of the pipeline using real data and shows with an example the application on clinical data for quantifying atrophy in Alzheimer's Disease using cortical thickness. We also compared the results with slices. The protocol of acquisition and demographics are fully detailed in (P. [START_REF] Bourgeat | AB burden in the temporal neocortex is related to hippocampal atrophy in non--demented elderly[END_REF]. The SC was used to register the surfaces to a template, which was obtained as an average of a subsample (20 individuals) randomly selected from the whole population.

Partially inflated surfaces

Using the 30 young healthy individuals randomly selected from the OASIS database [START_REF] Marcus | Open Access Series of Imaging Studies (OASIS): Cross--sectional MRI Data in Young, Middle Aged, Nondemented, and Demented Older Adults[END_REF], we first compared the results of the inflation step with Freesurfer. In this experiment, segmentations of the WM/GM interface were computed using Freesurfer. Thus, for each individual we obtained separated surface representations for left and right hemispheres. We evaluated the computation of PFS separately in terms of area preservation. A thorough comparison of different methods using local and global area, angle and length distortions was presented in [START_REF] Bonner | A quantitative comparison of three methods for inflating cortical meshes[END_REF]. In this paper we only report the measured local metric distortion in the 1--neighbourhood surrounding each mesh vertex. The local area distortion was defined for each vertex i as

∑ ∑ ∈ ∈ = T t T t A t A t A i M ~0 ~) ( ) ( ) ( ˆ (7)
and the length distortion as

∑ ∑ ∈ ∈ = 1 1 ~0 ) ( ˆN j ij N j ij D d d i M (8)
whereT ~is the set of triangles containing vertex i , 1 Ñ is the set vertices in the 1--

neighbourhood of i , ) (t A is the area of triangle t and ) ( 0 t A
is the area of triangle t on the original mesh. This is a measure of the average absolute change in area per triangle on the mesh. 0 in d and in d are the distances between vertex i and n on the original and modified surfaces, respectively. In our experiments, the implemented method for computation of PFS has been proven [START_REF] Bonner | A quantitative comparison of three methods for inflating cortical meshes[END_REF] Table 1 shows the average portion of vertices within this range for each method, computed for all 30 subjects in the sample space per lobe. Overall, the distance metric was distorted nearly equally by both CARET and FreeSurfer (23.03 ± 0.86% for CARET vs 23.11 ± 1.05% for FreeSurfer) whereas the area was distorted by 32.7 ± 1.42% using CARET vs 39.09 ± 1.71 with Freesurfer.

Shape context--based non rigid registration

The overall method was used to register the obtained surfaces to a common template.

For the experiments described using OASIS the template was obtained as a surface representation of the Colin atlas (D.L. [START_REF] Collins | Design and Construction of a Realistic Digital Brain Phantom[END_REF]. This template was obtained with

Freesurfer applying the segmentation pipeline to the initial Colin's T1--W MRI. For the experiments using clinical data a different template was obtained as an average surface representation of the population being studied.

In the first experiment using the OASIS database, we investigated the effect of the number of control points ( u ) employed to represent the simplified surfaces s M and s F on the registration quality. Mean absolute (MAD) and Hausdorff (HD) distances [START_REF] Gerig | Valmet: A New Validation Tool for Assessing and Improving 3D Object Segmentation[END_REF] between moving and target surfaces were used to assess the accuracy of our method to register two PFSs. Table 2 presents a summary of these results. By Increasing u from 500 to 1000 control points the accuracy accuracy in terms of mean absolute distance (MAD) was improved by 33% ( p <0.0001) and by 32% ( p <0.0001) when varying from 1000 to 2000 points for the left hemisphere and 23% ( p <0.0001) and 27% ( p <0.01) respectively for the right hemisphere. As expected, surface alignment is better when using more control points.

A trade--off between computational requirements and accuracy exists and must be defined depending on the application.

Overlap with Freesurfer labeling

In this experiment, we compared the localization of Freesurfer labeled regions (obtained as described in [START_REF] Desikan | An automated labeling system for subdividing the human cerebral cortex on MRI scans into gyral based regions of interest[END_REF])) after registration to the same template (Fig. 6). The Jaccard coefficient measured as the relation between the intersection and union of areas of labelled regions was used as the similarity metric. Overall, our method is comparable with the registration provided with Freesurfer. Averaging for the 30 individuals, in 6 regions of the left hemisphere and 6 of the right hemisphere the overlap (Jaccard) obtained with our method was higher than with Freesurfer. In other regions the accuracy was lower than with Freesurfer but in most of them these differences they were not statistically significant, except for the entorhinal cortex, middle temporal, pars triangularis and transverse temporal [START_REF] Desikan | An automated labeling system for subdividing the human cerebral cortex on MRI scans into gyral based regions of interest[END_REF]) in both hemispheres. Table 3 shows some examples of the jaccard coefficient averaged on the 30 individuals. The worst values were obtained in the smallest regions, where the overlap is highly affected by misregistration errors. Conversely in larger areas the overlap was fairly good, and the main folds were well aligned. It is important to notice that the statistical analysis will be performed later at a voxel basis. Fig. 6 illustrates an example of propagated labels between an individual and the template. It can be observed that although there is a high interindividual variability, the regions are well transferred through all the foldings representing different anatomical regions.

[Figure 6]

Clinical data: AIBL study

Following the procedure detailed in Section 2, consensus segmentations were obtained as described in [START_REF] Bourgeat | AB burden in the temporal neocortex is related to hippocampal atrophy in non--demented elderly[END_REF]) using nine atlases for initialization. Subsequently, cortical thickness maps were computed using our voxel--based approach [START_REF] Acosta | Automated voxel--based 3D cortical thickness measurement in a combined Lagrangian--Eulerian PDE approach using partial volume maps[END_REF]). The individuals' Automated anatomical labeling (AAL) (Tzourio--Mazoyer N, et al 2002) and cortical thickness were propagated to the resulting mesh and then to the average template after the SC registration.

As an example to describe the matching of similar regions of the brain, we randomly selected 26 individuals. Their consensus AAL labelings, obtained in the previous step, were propagated towards the template applying the whole pipeline. The labelings were then combined in a voting scheme to obtain a final AAL labeling in the template space (Fig. 9). After the registration, consistency in the coincidence of the major folding patterns was obtained across the individuals.

Cortical thickness comparison with Freesurfer

Cortical thickness maps from 50 individuals obtained with FS and milxCTE (Acosta et al. 2009) were compared. After registration and mapping to the common template, the thickness maps were regionally compared. Cortical thickness was comparable in all the regions. Similar patterns for thickness in gyri and sulci were obtained. Overall, a strong correlation between the results using both methods FS and milxCTE was found. In average, for all the regions the absolute computed difference was 0.23mm ( r =0.63) for the left hemisphere and 0.17mm ( r =0.65) for the right. Table 4 shows the results for some of the regions, left and right hemispheres.

[Figure 7]

[Figure 8]

Figure 7 shows the results of the correlations ( r ) for all the regions. Similar results were obtained in both hemispheres. Some examples of the individuals' thickness computed with milxCTE and FS are depicted in Fig. 8. As suggested in [START_REF] Klauschen | Evaluation of Automated Brain MR Image Segmentation and Volumetry Methods[END_REF]) Fresurfer tends to underestimate gray matter, therefore the computed thickness with our method tends to be higher. Although the slope on Fig. 8 is leaning towards the higher values of CTE, the correlations are still high as shown in table 4. In milxCTE the thickness estimation is entirely voxel--based with a sub--voxel initialization using partial volume maps (Acosta et al. 2009). This makes it very fast and accurate whereas FS is based on meshes deformations which may be computationally expensive. Not surprisingly, larger regions: Supramarginal (2400 vertices),

Inferior parietal (4200), superior temporal (2600), superior frontal (4500) presented the highest correlations as opposed to smaller regions such as rostral anterior cingulate (400), frontal pole (60) or caudal anterior cingulate (813) where small overlap differences may affect the measure.

Study of cortical thickness in Alzheimer's disease

We also investigated the ability of our method to detect cortical thickness differences between 81 healthy elderly individuals HC, and 32 AD patients. Fig. 9 shows an example of cortical thickness maps from two individuals, one from each group, propagated to the common template.

[Figure 9]

[Figure 10]

After applying the whole pipeline to all the individuals, vertex--wise t--tests were performed between NC--AD to identify regions where a significant atrophy existed. The cortical thickness values were corrected for age and False Discovery Rate (FDR) corrected p--values were obtained (5% threshold). In our experiments, the analysis was performed using the thickness maps with a Laplacian smoothing (5mm Full Width at half of the Maximum--FWHM) over the scalar cortical thickness maps to reduce the effects of discontinuities. The obtained differences between the two groups are illustrated in Fig. 10 and demonstrated qualitatively the consistency of the results. The results show significant differences between the two groups, with lowest thickness for the AD group with both methods. Significant differences were found in the hippocampus, parahippocampus, cingulate and in the temporal and frontal lobes between the two groups. In average, the absolute differences in cortical thickness for some regions were higher by using our method, as it can be seen in Figure 10 and in Table 4. These results corroborate what has been previously reported in the literature for Alzheimer's disease (Lerch and A. C Evans 2005), [START_REF] Querbes | Early diagnosis of Alzheimer's disease using cortical thickness: impact of cognitive reserve[END_REF].

Conclusion

In this paper we presented a new cortical surface processing pipeline, milxCTE, for accurate statistical regional analysis. It offers an alternative to traditional methods using surfac representations of the brain such as Freesurfer. It involved several steps, namely MRI segmentation, surface generation, topology correction, inflation and surface matching, yielding a representation of a whole population in a common space. In addition to the shape context which provides a global description at a local level, we included additional local landmarks, such as the sulcal depth map, to iteratively adjust the registration, yielding a meaningful matching between lobes. An update of the binary distribution of milxCTE and milxView is available 10 .

We demonstrated the validity of the method for use in clinical studies, by evaluating each step separately on real data, and then comparing the overall technique against Freesurfer.

This comparison showed that the labeling of major folding patterns is preserved and the cortical thickness computed by using our method is regionally consistent. In average, with

Freesurfer the obtained values of cortical thickness were lower than with milxCTE for single individuals. This is not a surprising finding as it has been previously shown and reported in larger studies. Our method performs the computation of cortical thickness in the voxel domain, which is very fast. Additionally, it uses advantageously the partial volume information to initialize at a subvoxel level, thereby yielding a very accurate result. However, a good correlation between both measures were found and different experiments comparing averaged cortical thickness between two populations allowed to measure the relative differences of atrophy. This preliminary study on clinical data showed regional differences between healthy elderly individuals and Alzheimer's disease patients. The most significant atrophy was measured in the temporal lobe, which is consistent with the published literature.

We intend to perform longitudinal clinical studies on Alzheimer's disease and other neurological disorders. As similar results were found in regional comparisons with Freesurfer, our method represents an alternative for the study of cortical thickness estimation across a population.

As the method allows the mapping of a whole population in a common space, in future topology is corrected on each separated hemisphere (5), partial inflation is performed (6) and eventually surfaces are registered towards a common coordinate system (7).

  Fig. 4 (Top) shows the different steps in the evolution of the PFS with the sulcal depth maps computed at each step. The color coded sulcal depth map allows the visual localization of the main folding patterns. The sulcal depth map is computed as the cumulative path length at each point i as

  Fig 4 Bottom depicts the comparison of Area distortion after inflation using CARET and Freesurfer.

  the best trade--off for a good matching of the folds. Given the individual costs ] of points, the next step is to find the perfect matching by minimizing the total cost of the bijective correspondences

Freesurfer

  (A.M. Dale, B. Fischl, and M.I. Sereno 1999), (B. Fischl, M. I Sereno, and A. M Dale 1999), (B.Fischl and A. M Dale 2000). All the methods were implemented in C++, incorporated in a plugin called milxCTE devised to compute the cortical thickness as described in(Acosta et al. 2009), and is part of our downloadable software platform milxView 5 , and TAGLUT (Topological And Geometrical Library --a Useful Toolkit) 6[START_REF] Favreau | TAGLUT: Topological And Geometrical Library --a Useful Toolkit[END_REF] , which provides the topological and geometrical tools for genus zero mesh generation. Milxview utilises the open source ITK 7 and VTK 8 libraries. Access Series of Imaging Studies (OASIS) database[START_REF] Marcus | Open Access Series of Imaging Studies (OASIS): Cross--sectional MRI Data in Young, Middle Aged, Nondemented, and Demented Older Adults[END_REF]) 9 , we randomly selected 30 young healthy individuals. The scans were T1--weighted Magnetization Prepared RApid Gradient Echo (MP--RAGE) in sagittal orientation with isotropic 1mm 3 resolution (256x256x128 pixels). For the first set of experiments, inflation and registration parts, the segmentations of the WM/GM interface were computed using Freesurfer. Thus, for each individual we obtained separated surface representations for left and right hemispheres. The SC was used to register the surfaces to a common template, which was obtained after segmentation of the Colin atlas(D.L.[START_REF] Collins | Design and Construction of a Realistic Digital Brain Phantom[END_REF] with Freesurfer. The resulting meshes contained an average of 300000 vertices. 5 http://aehrc.com/biomedical_imaging/milx.html3.1.2 Clinical data in Alzheimer's disease: the AIBL studyFor the application of the proposed method to the study of Alzheimer's disease, MR images from AD patients, individuals with mild cognitive impairment (MCI), and HC were included. They were randomly selected from the individuals enrolled in the Alzheimer Imaging Biomarkers and Lifestyle longitudinal (AIBL) study (Kathryn A[START_REF] Ellis | The Australian Imaging, Biomarkers and Lifestyle Study of ageing (AIBL): methodology and baseline characteristics of 1112 individuals recruited for a longitudinal study of Alzheimer's disease[END_REF]). The AD participants met NINCDS--ADRDA criteria for probable AD[START_REF] Mckhann | Clinical Diagnosis of Alzheimer's Disease: report of the NINCDS--ADRDA Work Group under the auspices of Department of Health and Human Services Task Force on Alzheimer's Disease[END_REF] while all the MCI subjects met the Petersen criteria[START_REF] Petersen | Mild cognitive impairment: clinical characterization and outcome[END_REF] of subjective and objective cognitive difficulties, in the absence of dementia or significant functional loss. The remaining participants in the study were healthy elderly volunteers. All patients were recruited from the Austin Health Memory Disorders and Neurobehavioural Clinics. All of the subjects were scanned with a standardized protocol. Thus, sagittal T1 weighted MR images were acquired using a standard 3D MPRAGE sequence at 3T, with in--plane resolution 1x1mm, slice thickness 1.2mm, TR/TE/T1=2300/2.98/900, flip angle 9° and field of view of 240x256 voxels and 160

  to perform better than FreeSurfer (B. Fischl, M. I Sereno, and A. M Dale 1999). Fig 4 Bottom depicts a comparison of the local area distortion after inflation using CARET and Freesurfer. Averaging over the 30 individual brains surfaces, it produced substantially fewer area and distance distortions. 45.89% of vertex no distortion), compared to 34.83% for FreeSurfer.

Figure 1 :

 1 Figure 1: Overall pipeline for cortical mapping with milxCTE. After WM/GM segmentation (1--4),

Figure 2 :

 2 Figure 2: Cutting and patching of non--separating loops. a) Computation of the cut locus B associated to a point p on a 1--genus surface M. b) Cut locus B associated to a point p on a 2-genus surface M. c) Computation of the shortest non--separating loop associated to a point on the path S B , subpart of the cut locus B associated to the point p.

Figure 3 :

 3 Figure 3: Example of topological correction of a,b) a tunnel, c,d) two handles.

Figure 4 :

 4 Figure 4: Top: Original cortical mesh is flattened to a PFS and color coded sulcal depth. (Eq. 2) is mapped over each point in the surface as inflation progresses (20,50 and 100 iterations). a) Original Mesh b) 20 iterations c) 50 iterations d) 100 iterations. Bottom : Local Area distortion map after inflation using CARET and Freesurfer.

Figure 5 :

 5 Figure 5: Finding correspondences between two individuals. After PFS are obtained, the surface M is registered towards the fixed F target. Ms and Fs are the control points to match.

Figure 6 :

 6 Figure 6: Propagated labels from the template d) to the patient c) after the registration of both PFS surfaces a), b).

Figure 7 :

 7 Figure 7: Regional mapping of the correlation ( r ) between cortical thickness computed with both methods: milxCTE and Freesurfer, left and right hemispheres.

Figure 8 :Figure 9 :

 89 Figure 8: Example of correlations in some regions left/right hemispheres. a,b) SFG: superior frontal gyrus and c,d) IPG: Inferior parietal gyrus [SFG (left H), r =0.80 ] [SFG (right H), r =0.75] [IPG (left H), r =0.86] [IPG (right H), r=0.86] (corrected values)

Figure 10 :

 10 Figure 10: Top: Vortex--wise differences in averaged cortical thickness between the two groups AD and HC. Bottom: Statistical p--values Map after two sampled t--test.

Figure 1 Figure

 1 Figure 1

Table 1 :

 1 TABLES Average percentage of vertices with A M ˆ in range [0.8,1.2]

	Hemisphere		CARET	FreeSurfer
	Left	Frontal	49.31	38.67
		Temporal	45.77	37.69
		Occipital	39.94	31.92
		Parietal	46.67	32.86
		Total	42.41	32.94
	Right	Frontal	49.59	39.00
		Temporal	45.13	37.49
		Occipital	39.50	32.91
		Parietal	46.76	32.44
		Total	42.47	32.91
	Table 2: Averaged Hausdorff Distance (HD), Mean Absolute Distance (MAD) and area
	relation (% A M --F ) betweeen moving M and fixed F surfaces after shape context
	registration. Mean (standard deviation), as the number of sampling control points varies. 30 left
			and 30 right hemispheres.
			left h			right h
	NPts	500	1000	2000	500	1000	2000
	HD [mm]	8.45 (1.9) 7.25 (1.97) 6.29 (1.71)	14.97	13.24	13.11
					(7.13)	(7.83)	(8.25)
	MAD [mm] 0.74 (0.05) 0.49 (0.04) 0.35 (0.03) 1.17 (0.58) 0.86 (0.61) 0.74 (0.83)
	% A	98.08%	98.45%	98.72%	97.02%	97.72%	98.8%

Table 3 :

 3 Averaged (standard deviation) jaccard coefficient on some anatomical labels computed after registration with Shape Context (SC) and with Freesurfer (FS).

	SC	FS

Table 4 :

 4 Average (standard deviation) cortical thickness, computed with milxCTE and

	Freesurfer and correlation between both methods in the AAL Regions.
		milxCTE		Freesurfer	Comparison
	Structure		St Dev		St Dev		r
		Average(mm)		Average(mm)		Difference
						(mm)
	Supramarginal L	2.38	0.20	2.45	0.13	0.07	0.69
	Supramarginal R	2.38	0.17	2.45	0.14	0.07	0.65
	Superiorfrontal L	2.49	0.17	2.59	0.12	0.01	0.64
	Superiorfrontal R	2.54	0.20	2.54	0.11	0	0.73
	Rostralmiddlefrontal L	2.32	0.18	2.28	0.11	--0.04	0.72
	Rostralmiddlefrontal R	2.36	0.18	2.25	0.10	--0.11	0.67
	Rostralanteriorcingulate	2.62	0.32	2.64	0.21	0.02	0.68
	L
	Rostralanteriorcingulate	2.67	0.29	2.52	0.18	--0.15	0.78
	R
	Precuneus L	2.53	0.27	2.28	0.15	--0.25	0.79
	Precuneus R	2.45	0.26	2.28	0.14	--0.17	0.8
	Precentral L	2.21	0.28	2.43	0.13	0.22	0.72
	Precentral R	2.26	0.27	2.38	0.15	0.12	0.83
	Posteriorcingulate L	2.57	0.24	2.38	0.14	--0.19	0.57
	Posteriorcingulate R	2.49	0.21	2.33	0.12	--0.16	0.72
	Pericalcarine L	2.29	0.29	1.74	0.12	--0.55	0.65
	Pericalcarine R	2.2	0.25	1.76	0.12	--0.44	0.72
	Pars triangularis L	2.41	0.15	2.32	0.11	--0.09	0.65
	Pars triangularis R	2.3	0.16	2.34	0.12	--0.04	0.62
	Inferior temporal L	2.87	0.36	2.67	0.15	--0.2	0.6

http://www.fil.ion.ucl.ac.uk/spm/

http://surfer.nmr.mgh.harvard.edu/

http://research.ict.csiro.au/software/milxview

http://www.cma.mgh.harvard.edu/ibsr

http://www.jmfavreau.info/?q=en/taglut

http://www.itk.org/

http://www.vtk.org/

www.oasis-brains.org

http://aehrc.com/biomedical_imaging/milx.html

work we will perform multi--feature comparison of additional imaging biomarkers for Alzheimer's disease. Thus, not only cortical thickness is compared but also other features from different modalities such as sulcal depth, WM/GM MRI contrast, GM volume, β --Amyloid burden or WM integrity in the same vertex--wise common space.

Algorithms


Jérémy Lasbleiz Ab

Hervé Saint-Jalmes

Régis Duvauferrier

Anita Burgun

Creating a Magnetic resonance imaging ontology

Keywords: MRI, MRI Simulator, OWL, ontology

The goal of this work is to build an ontology of Magnetic Resonance Imaging. The MRI domain has been analysed regarding MRI simulators and the DICOM standard. Tow MRI simulators have been analysed: JEMRIS, which is developed in XML and C++, has a hierarchical organisation and SIMRI, which is developed in C, has a good representation of MRI physical processes. To build the ontology we have used Protégé 4, owl2 that allows quantitative representations. The ontology has been validated by a reasoner (Fact++) and by a good representation of DICOM headers and of MRI processes. The MRI ontology would improved MRI simulators and eased semantic interoperability.

Introduction

Magnetic Resonance Imaging is the most versatile diagnostic imaging technique. It can study T1, T2, diffusion, PH, temperature, spectroscopy… of tissues and of course make images. The vocabulary used by medical imaging constructors is very heterogeneous [START_REF] Gibaud | The quest for standards in medical imaging[END_REF] and physical phenomena involved during MRI are very complex. So the MRI domain needs ontology to make the MRI community sharing the same concepts. To build our ontology we will take into account two MRI representations: MRI simulators and DICOM. The DICOM is an applicative representation with daily-use concepts. MRI simulators give a representation of complex physical phenomena that are involved in MRI and that are not describe in DICOM. The fusion of MRI simulators and DICOM concepts is needed to represent MRI examinations not only in an administrative way but in a useful way for radiologist interpretations.

Material and methods

Analyzing DICOM [2]

The DICOM standard is divided in different parts. The relevant part for MRI is C.8.13 « Enhanced MR Image ». It is a section of the standard part 3: « Information Object Definition ». All concepts of this part, and their DICOM tags, will be included in our ontology, thus will give a semantic interoperability to the ontology. In DICOM, there is a lack of definition for an ontology. We will fill this gap by domain expert definitions, thanks to MRI simulators analysis.

Analyzing MRI simulators

We decided to analyze two MRI simulators JEMRIS and SIMRI. SIMRI [START_REF] Benoit-Cattin | The SIMRI project : a versatile and interactive MRI simulator[END_REF], is implemented in C language and is based on the Bloch equations. It enables simulations of 1D, 2D, and 3D images. Although simple, the user interface requires the use of C. The simulator is divided in different parts: -Model (Virtual object): Each voxel of the virtual object contains a set of physical values that are necessary to compute the local spin magnetization vector with the Bloch equations. These values are the proton density and the two relaxation constants T1 and T2.

-MRI sequence: During an MRI experiment, the object is placed in a static magnetic field B0 and is excited by electromagnetic events of two types: RF pulses (B1 field) and magnetic field gradients.

-The acquisition of the object magnetization state is stored as a complex signal in the k-space to obtain the image. This part is divided in 4 parts: The free precession, precession with application of gradients (specified by its duration and the gradient magnitudes in the three spatial directions), signal acquisition (number of points to capture, bandwidth, readout gradient magnitude and position of this signal in the kspace), the application of RF pulses (specified by its duration, a flip angle and the rotation axis). RF inhomogeneity and gradient non-linearity are not simulated. The user can define the echo train and sequence parameters (repetition time, echo time, flip angle…). Chemical shift and susceptibility artefacts are modeled.

JEMRIS [4-5]

, is a C++ software with XML tags. It uses an optimized library for numerical solutions equations needed to simulate complex RF-pulses. It can deal with multichannel Tx-Rx coil geometries and configurations, nonlinear gradients, chemical shift, reversible spin dephasing (T2*), susceptibility-induced off-resonance, temporal varying processes of the object (e.g., movement or flow), and concomitant gradient fields. The graphical user interface (GUI) is divided in three: one for interactively designing the MRI sequence, another for defining the coil configuration, and one for the setup and execution of the main simulator. The software is divided in 5 classes: sample (describes the physical properties of the object) signal (holds information about the MR signal) model (describes the functionality for solving the physical problem) coil (contains the code for spatially varying RF transmission and signal reception), sequences. The sequence loop is represented as a left-right ordered tree with loops (Fig1). The xml language has been used to serialized C++ objects, describing the different steps of each sequence. The management of time interval has also been taken into account and formalised. The different modules interact with each other.

Fig.1: Echo Planar Imaging sequence schema in JEMRIS [4] yellow =loops , blue =

pulses, green = intervals.

Using Protégé 4, owl2, Ontology validation

To build our ontology, we will use Protégé [6], which is a free, open source ontology editor and knowledgebase framework and the owl language. In our case, the domain has a lot of quantitative informations so we choose owl2, which allows us to define quantitative data properties. First of all, we have taken into account concepts from DICOM and secondly we have added concepts from MRI simulators. We will use an ontology classifier FACT++ to check the ontology consistency. The ontology will be validated by the analysis of 10 MRI examination DICOM headers, extracted with OSIRIX [7] and the possibility to write sequences with the ontology.

Results

Ontology taxonomy

The main classes of ontology taxonomy are: Object of the study: Defined by its size, voxels size, properties (T1, T2, Proton Density, Diffusion, Contrast enhancement cinetic), T2*, movements (general and flux) Device: Magnet (intensity, shape, kind) coil (receiver coil, transmitter coil, multielement coil, region) Gradient (magnetic field, slice selection, diffusion…) Sequence, from this point our vision is different from JEMRIS. Actually, the representation of loop in a vertical way (fig. 1) of physical events that are horizontal (dependant to time) and independent cannot be included in an ontology. Therefore we have divided sequences in elementary events: radiofrequency pulse, slice selection gradient, readout gradient… according to SIMRI description of events. The signal acquisition modeled, has to be formalised by a mathematical way thanks to Bloch equations resolution as in the two softwares. The formula will be integrated in the ontology. Acquisition results will be divided in: image, quantitative result… Organisation of sequences in taxonomy is a difficult management. An article [START_REF] Ge | An Interactive Taxonomy of MR Imaging Sequences[END_REF], written in a didactical goal, has organized sequences with their technical characteristics and with loops. Taxonomy doesn't have loop and the problem is that sequences can be a mix of different techniques that can't be organised in taxonomy. The solution we have chosen is to classify sequence according to their goals. This solution is intuitive for clear goal: diffusion, angiography image… but less obvious for contrast sequences. So we have chosen to start with a general taxonomy of sequences (Fig. 2), adding to each of them the Weighting of final images: T1Weighted, T2Weighted, DPWeighted and T2*Weighted.

Constructor acronyms of sequences have been added as synonyms of sequence name. We extract DICOM headers of 10 MRI examinations with OSIRIX Métadonnées. The analysis shows that concepts of DICOM headers are well represented in the ontology. The problem is that MRI constructors don't share the same DICOM tags for the same concept.

Discussion

There is only one work about MRI and ontology. It concerns brain functional MRI [START_REF] Nakai | Ontology for FMRI as a biomedical informatics method[END_REF] and are interested in all the process and not only MRI. However it has already shown the need of ontology in the domain. JEMRIS have also, by using XML, shown the interest of web semantic in physical process description. DICOM also need to be improved with definitions and rules that ontology could define. Our ontology can increase the semantic interoperability in MRI. An ontology has already be implemented on a PACS in that goal [START_REF] Dl. Rubin | Medical Imaging on the Semantic Web: Annotation and Image Markup[END_REF] but not for MRI examinations.

Fig. 2 :

 2 Fig.2: Contrast sequence Taxonomy Acquisition parameters are divided in two essential parts: parameters modifying image geometry and parameters modifying image contrast. The ontology relations are: Different kinds of relations between concepts will be defined: General: Has_a , Has_Parameters… ; Quantitative : Has_Value, Has_Unit… Owl2 permits quantitative representation of classes. The relations between classes are then: A Has_Modifyer B, A Increase_When_Decrease B, A Decrease_When_ Increase B will permit to describe variations of parameters.

Corresponding Author.