PREDICTING CHANGES IN VISUAL APPEARANCE OF PERIODIC SURFACE FROM BRDF… PREDICTING CHANGES IN VISUAL APPEARANCE OF PERIODIC SURFACE FROM BRDF MEASUREMENTS

Keywords: Surface patterning, glossy behaviour, high resolution BRDF measurements, diffraction effect

The present study focuses on the optical properties of functionalized surfaces and how the surface geometry impacts them. Physical measurements of reflected light are required to understand the visual aspect of such surfaces. Bidirectional Reflection Distribution Function (BRDF) is evaluated in order to identify and understand physical effects induced by different surface functionalization. BRDF measurements of high angular resolution allow us to observe diffraction phenomenon at an unusual scale. Experimental results were compared to theoretical calculation, and then diffraction effect was confirmed. This study stresses that diffractive phenomenon could have an impact on visual aspect even for simple pattern geometry. In the same way, if more complex periodic patterns are considered, as multistate periodicity pattern for example, it could rise to important modification of surface visual aspect.

Introduction

Carefully chosen surface modification can provide new properties for applications in wetting, haptics, optics to almost any surface [START_REF] Xia | Bio-Inspired, Smart, Multiscale Interfacial Materials[END_REF]. Yet, surface modification can also impact visual properties and thus aesthetics of the product. Visual appearance is a major criterion for customers: the final aspect of surface could be an issue for surface patterning. Link between surface pattern, physical measurement of light reflection such as BRDF (Bidirectional Reflection Distribution Function) and visual appearance is not trivial [START_REF] Flys | Characterization of surface topography of a newly developed metrological gloss scale[END_REF]. The goal of our study is to bring first pieces of knowledge in the comprehension of the relation between surface patterning, BRDF and visual appearance of such functionalized surfaces. In order to understand all phenomena at stake, high angular resolution BRDF measurements are required to be compared to theoretical calculations based on grating formula.

Samples description

In the sake of simplicity, optical properties of hydrophobic surfaces are studied on the example of model hydrophobic patterns with simple geometry (Fig. 1). Each sample consists of a glass slide with the patterned hybrid silica layer on its top (Fig. 1 a). Used patterns were cylindrical micro pillars (height (h) and diameter (d): 10 µm) organized within hexagon lattice. The lattice constant (a) varied from 20 to 80 µm (Fig. 1b), which in terms of the surface coverage rate (surface parameter defining the wetting properties [START_REF] Gauthier | Role of Kinks in the Dynamics of Contact Lines Receding on Superhydrophobic Surfaces[END_REF]) represents the variation within the range between 1.4% (80 µm lattice) and 22.7% (20 µm lattice). In case of hexagonal lattice this latter parameter is calculated by the following formula:

𝜌 = 𝜋 𝑑² 𝑎² √3 6 (1)
Residual sol-gel layer between the glass surface and pillars is necessary to provide the adhesion of pillars to the surface. 300 nm thickness of this layer was chosen to prevent from optical guiding effects.

Device fabrication process consists of two main steps: silica sol -gel chemistry and nanoimprint (for detailed description see [START_REF] Dubov | Superhydrophobic silica surfaces: fabrication and stability[END_REF]). A hybrid silica sol-gel is prepared and is spin coated on a glass substrate (5 cm x 5 cm in our case). Then the nanoimprint step is applied: a PDMS (Polydimethylsiloxane) mold, previously patterned, is put on the coated glass substrate. The pattern is then replicated from the mold to the sol -gel layer under controlled conditions (1 bar pressure and temperature of 110°C during 45 minutes). This fabrication process all ows creating various types of functionalized surfaces on glass substrate, with high reproducibility.

All samples were examined in a "Spectralight QC X-Rite" light booth, under D65 illumination. Surface patterning only slightly modifies transmission properties of glass-based samples. Nevertheless there are significant visual differences in reflection between samples with relatively high (>8%, a=20 and 30 µm) and lower surface coverage (5%, a40µm) and especially at grazing angles (Fig. 2). Meanwhile, samples with higher coverage rate (>8%, a=20 and 30 µm) look diffuse-like in reflection, as they exhibit a white-coloration typical for lambertian-like rough surfaces. Samples with low surface coverage (5%, Fig. 2a) give the impression of smooth and specular-like surfaces (as illustrated on Fig. 2b, there is no blurring on the image of coloured grid reflected by such surfaces) . These observations suggest that surface coverage rate of such patterned surfaces impact their gloss .

In order to understand visual rendering of samples, optical properties have to be examined. Integral entities such as gloss and haze are not sufficient and s patially resolved measurements should be applied to understand this effect. Due to the highly glossy behaviour of samples (especially for low surface coverage), a goniospectrophotometer with high angular resolution is required to perform relevant measurements. As light scattering phenomena is typically represented as a function of polar and azimuthal angles, an angular convention, used to describe illumination and detection angles on this paper, is presented on Fig. 3, adapted from [START_REF] Hetland | Numerical studies of the scattering of light from a two-dimensional randomly rough interface between two dielectric media[END_REF]. iandidescribe illumination direction whilerandrdescribe reflected direction.

We now turn to the light scattered from the samples introduced in the previous section. The phenomenon of light scattering from periodic surface structures is known as diffraction, and it has been studied for centuries, and it is still of contemporary interest in science and engineering [START_REF] Born | Principles of optics: electromagnetic theory of propagation, interference and diffraction of light[END_REF]. When incident light is scattered from a periodic struct ure, known as a grating, the reflected light typically consist of a set of pronounced peaks in the reflected intensity distribution. Each of these peaks is referred to as diffractive orders, and the positions of these peaks are mainly determined by the wavelength of the incident light, the type of grating and parameters defining it.

Mathematically these angular positions are determined by the grating formula [START_REF] Loewen | Diffraction gratings and application (25-53[END_REF] that we now will introduce (in a pedagogical fashion). To this end, we will assume the scattering geometry depicted in Fig. 3. The planar surface supporting the pillars is assumed to coincide with the x1x2 plane. In the region x3 < 0 we have a medium that is characterized by the dielectric constant , that is assumed to be real and positive (i.e. a dielectric); and in the region that is x3 > 0 and outside the pillars, we assume to have vacuum (0 = 1). Into this structure, light of wavelength  is incident at the angles of incidence (i, i). The component of the incident wave vector k that is parallel to the planar surface of the substrate is then

k || = 𝜔 𝑐 𝑠𝑖𝑛𝜃 𝑖 (𝑐𝑜𝑠ϕ 𝑖 , 𝑠𝑖𝑛ϕ 𝑖 , 0), (2a)
where we have introduced the frequency of the incident light, and c denotes the speed of light in vacuum [/𝑐=/]. Similarly, the lateral wave vectors of a beam reflected into the angles of reflection (r, r) are (see Fig. 3)

q || = 𝜔 𝑐 𝑠𝑖𝑛𝜃 𝑟 (𝑐𝑜𝑠ϕ 𝑟 , 𝑠𝑖𝑛ϕ 𝑟 , 0), (2b)
The lattice of the periodic structure is invariant under a translation by the vector x|| (l) =l1a1+l2a2, where ai (i = 1; 2) are the two non collinear primitive translation vectors (in the x1x2 plane) that characterize the lattice (see Fig. 1), l=(l1,l2) and li is an integer (li ∈ ℕ). The area of the unit cell is Ac=|a1 x a2|. Associated with the direct lattice, characterized by a1 and a2, is the reciprocal lattice defined by the translation vectors where m=(m1,m2), mi ∈ ℕ, and the primitive translation vectors, bi, are defined from the relation (i, j = 1, 2) 𝒂 𝑖 . 𝒃 𝑖 = 2𝜋𝛿 𝑖𝑗 , [START_REF] Dubov | Superhydrophobic silica surfaces: fabrication and stability[END_REF] where ij represents the Kronecker delta.

𝑮 || (𝒎) = 𝑚 1 𝒃 1 + 𝑚 2 𝒃 2 , (3)
In this work we assume that the direct lattice formed by the pillars is of the hexagonal type, so that the reciprocal lattice also is hexagonal. Without loss of generality, we orient the coordinate system so that the primitive lattice vectors may be taken as

𝒂 1 = 𝑎 √3 2 𝒙 ̂1 + 𝑎 1 2 𝒙 ̂2 𝒂 2 = -𝑎 √3 2 𝒙 ̂1 + 𝑎 1 2 𝒙 ̂2 , (5)
where a > 0 is the lattice constant and a caret over a vector indicates a unit vector. The corresponding primitive translation vectors for the reciprocal lattice, bi, are determined from Eq.( 4).

After the incident light, characterized by the wave vector k|| [Eq. ( 2a)], has interacted with the periodic structure, the diffracted wave of order m in either reflection or transmission will have the wave vectors which component is parallel to the surface of the substrate

𝒒 || (𝒎) = 𝒌 || + 𝑮 || (𝒎) . (6)
Equation ( 6) is the grating formula that we seek, and it works equally well in both reflection and transmission. It should be noticed that for a given diffractive order m, the wave vectors 𝒒 || (𝒎) does To relate an open diffractive channel in reflection to the angles of reflection (r (m) , r (m) ) that defines it, we have to equate the right-hand sides of Eq. (2b) and Eq. ( 5)

not
𝒒
In obtaining Eq. ( 7) we have used that 𝑞 || = ( 𝜔 𝑐 ⁄ )𝑠𝑖𝑛𝜃 𝑟 , and the two argument arctangent function, denoted arctan(y,x) is used in order to determine the appropriate quadrant of the resulting angle (and to properly treat the case when 𝑞 1 (𝒎) = 0). In most computer languages this function is named atan2(.,.).

Materials and methods

BRDF measurements were realized with ConDOR (Conoscopic Device for Optical Reflectometry), a goniospectrophotometer with very high angular resolution (0.018°) dedicated to BRDF measurements within the specular zone. Measurements were performed using spatial multiexposure and time multiexposure. Spatial multiexposure allows to enlarge the range of detection of ConDOR setup from 2° (±1° around the specular direction) to 8° (±4°). Time multiexposure method provides high dynamics range (due to dealing with saturation issues of the CCD detector) and thus enables to detect the entire specular peak. Complete details on ConDOR goniospectrophometer as well as spatial and time multiexposure techniques are published elsewhere [START_REF]SPIE 9018, Measuring, Modeling, and Reproducing Material Ap pearance. Evaluation of the shape of the specular peak for high glossy surfaces[END_REF].

Samples with following covering rates were measured: 22.7%, 10.1%, 3.6% and 1.4%. Angular configurations of illumination and detection were fixed to: i=r=45°,i=r=0°. A V() filter was used during the measurements in order to reproduce the spectral sensitivity of human visual perception.

Gratings formula, described in the previous section, was implemented in order to predict the angular position of the diffractive orders. Illumination parameters were fixed to i=45° and i=

Results

BRDF measurements with ConDOR

All BRDF measurements exhibit a common structure of the spatial distribution of reflected light as shown in Fig. 4: a central specular peak surrounded by several secondary peaks, organized according to a hexagonal lattice. Angular distance between specular peak and secondary pea ks varies as a function of the pattern's size: a large period pattern induces a small angular distance, as presented in Fig. 5. Plane cuts presented in Fig. 5 were realized by representing the plane (red line) defined by the position of the specular peak and of one secondary peak (yellow crosses). For a=20 µm and a=80 µm, angular distances between the specular peak and its nearest neighbor secondary peak are respectively, 2.1° and 0.56°. The maximum of the specular peaks for all BRDF measurements are reported on Fig. 6a. BRDF values vary from 5.91 sr -1 (for the densest pattern) to 241.6 sr -1 (the least dense pattern). As expected low surface covering rate leads to high level of BRDF. BRDF values of secondary peaks are also influenced by pillar covering rate. Fig. 6b represents the ratio between BRDF values of specular peak and the averaged secondary peaks as a function of covering rate. The increase in covering rate induces the decay in this ratio, meaning that secondary peaks become more and more intense (this can be also observed on Fig. 5). Pattern with s mall covering rate reflects light close to the specular direction (in a +/-1° window), and the intensity of the specular peak makes the major contribution to BRDF compared to secondary peaks. While in case of the pattern with high covering rate the light is diffused in a wider angular range and the intensity of secondary peaks increases significantly.

Theoretical calculation results and comparison to measurements

As mentioned in Section 2, spatial distribution of the light reflected from the samples contains different diffractive orders defined by two angles (polar and azimuthal angles). Those angular positions are shown in Fig. 7a and Fig. 7b, respectively for patterns with a=20 μm and 80 μm.

As expected, those calculations predict a hexagonal organization of diffractive orders, and a strong relation between the pattern period and the angular distribution of diffractive orders. For instance, angular distance between two neighbors diffractive orders for a=20 μm pattern and 80 μm pattern are respectively 2.16° and 0.544°.

Measurements and calculations are compared on the example of the sample with 3.6% covering rate (Fig. 8). Fig. 8a shows a good correlation between measurements (yellow dots) and calculations (blue crosses) and confirms the presence of diffraction phenomenon. However, a small misalignment is observed resulting in small angular shift between experimental data and calculations. It is corrected in Fig. 8b through the 5° turn in azimuthal angle of illumination (i=5°). We attribute it to the inevitable misalignment of the micrometric pattern of the sample with the respect to the incidence plane. We finally end up with a very good angular correlation presented on Fig. 8b.

Discussion

Large gratings (a=80 μm in our case) induce an intense specular peak (241 sr -1 ), surrounded by several secondary diffractive peaks with a weaker intensity (2500 times smaller). All those peaks are localized in a ±0.75° window around the specular direction. As it is schematically shown in Fig. 9a, such patterns should enhance gloss property, due to their highly localized response in the specular direction. Dense gratings ( a=20 μm for this study) show a different behaviour. A specular peak (5.9 sr -1 ) is surrounded by only 40 times smaller secondary peaks which appear within a larger (±4°) window. With the same schematic representation, shown in Fig. 9b, the gloss zone is both enlarged and attenuated. Thereby, such diffractive pattern enhances the non-specular reflection of light resulting in the enhanced haze property. Finally, the variation of angular separation between specular peak and secondary peaks (as well as the variation in the relative intensity of diffractive orders) with the lattice period raise questions on the definition of gloss zone.

The diffraction phenomenon is rarely considered in the analysis of visual appearance. However, the experimental data presented in this paper prove that diffraction effect should not be neglected in case of surface functionali zation with periodic patterns (rather common for nanoimprint patterning technique). Moreover, patterns described here present a simple organization and the more complex lattice shape or multi -scale patterns could induce a stronger optical response and may significantly modify surface visual aspect. Furthermore, as rendering models based on microfacets theory do not account for light diffraction, they will fail to produce realistic visual rendering at least for the surfaces with periodic patterns.

Conclusion

Functionalization through surface patterning may significantly influence the visual properties of the object, especially in case of transparent objects. In the case of surface patterning with periodic arrays of micropillars, spatially resolved measurements showed the presence of numerous diffraction orders with angular separation and relative intensity stretching with the increase of surface density of pillars. Diffractive effects can induce slight modifications of visual aspect of such surfaces. Large pattern period tends to enhance gloss property of surface due to the localized angular distribution of diffractive peaks in the specular direction; while pattern with smaller period induces a larger angular distribution leading to a more diffuse surface aspect. Furthermore, those effects could be dramatically enhanced when dealing with more complex lattice patterns, consequently diffraction phenomenon must not be neglected when visual aspect of periodic system is studied.

Figure 1 -

 1 Figure 1 -Schematic description of the geometry of samples: a) stack of layers (side view), b) top view on hexagonal gratings of different lattice constant (20 and 80 µm)

Figure 2 -

 2 Figure 2 -Photographs of examined samples taken in a light booth (Spectralight QC X-Rite) under D65 illumination at grazing angles: a) comparison of samples with the highest (22.7%) and the lowest (1.4%) available surface covering rates; b) comparison of sample with the highest surface covering rate with samples with intermediate surface covering values.

Figure 3 -

 3 Figure 3 -Geometrical description of illumination and reflected angles: iandidescribe illumination direction andrandrdescribe detection direction.

  necessarily describe a propagating mode (an "open channel"); only those modes in reflection that satisfies 𝑞 || (𝒎) = |𝒒 || (𝒎) | < 𝜔 𝑐 ⁄ are propagating (in the ambient medium), and only when 𝑞 || (𝒎) = |𝒒 || (𝒎) | < √𝜀𝜔 𝑐 ⁄ is a mode propagating in transmission (assuming  positive).

  || = 𝒒 ||(𝒎) and solve the resulting equation system for those angles. The results are

Figure 5 -

 5 Figure 5 -Comparison of BRDF data from samples with different lattice constants (a=20 and 80 µm). Plane cuts are given in logarithmic scale.

Figure 4 -

 4 Figure 4 -BRDF measurement of sample a=50 µm at logarithmic scale: a) top view; b) a 3D representation. Experiment performed at i=r=45°,i=r=0° under unpolarized illumination with V() filter.

Figure 7 -Figure 6 -

 76 Figure 7 -Calculated angular positions of diffractive orders for a=20 μm (a) 80 μm (b). Calculation performed with parameters of the experiment: i=45°,i=0° and =555 nm

Figure 8 -

 8 Figure 8 -Positions of calculated (blue crosses) and experimentally observed (yellow points) diffractive orders: a) Illumination configuration are i =45° and i=0 for both measurements and calculations; b) calculation at i =45° and i=5° (to reproduce the misalignment of the pattern).

Figure 9 -

 9 Figure 9 -Schematic representation of the gloss behaviour of patterned samples with different periods: a) large period and b) small period.

Acknowledgement

Authors express deep gratitude to A.M. Rabal (LCM, LNE/CNAM, Trappes, France) for her precious help with ConDOR measurements.


Sylvain Takerkart ⋆a
email: sylvain.takerkart@univ-amu.fr

Guillaume Auzias

Lucile Brun

Olivier Coulon

Structural Graph-Based Morphometry: a multiscale searchlight framework based on sulcal pits

Keywords: morphometry, brain, sulcal pits, graph kernel, searchlight, multi-scale methods

Studying the topography of the cortex has proved valuable in order to characterize populations of subjects. In particular, the recent interest towards the deepest parts of the cortical sulci -the so-called sulcal pits -has opened new avenues in that regard. In this paper, we introduce the first fully automatic brain morphometry method based on the study of the spatial organization of sulcal pits -Structural Graph-Based Morphometry (SGBM). Our framework uses attributed graphs to model local patterns of sulcal pits, and further relies on three original contributions. First, a graph kernel is defined to provide a new similarity measure between pit-graphs, with few parameters that can be efficiently estimated from the data. Secondly, we present the first searchlight scheme dedicated to brain morphometry, yielding dense information maps covering the full cortical surface. Finally, a multi-scale inference strategy is designed to jointly analyze the searchlight information maps obtained at different spatial scales. We demonstrate the effectiveness of our framework by studying gender differences and cortical asymmetries: we show that SGBM can both localize informative regions and estimate their preferred spatial scales, while providing results which are consistent with the literature. Thanks to the modular design of our kernel and the vast array of available kernel methods, SGBM can easily be extended to include a more detailed description of the sulcal patterns and solve different statistical problems. Therefore, we suggest that our SGBM framework should be useful for both reaching a better understanding of the normal brain and defining imaging biomarkers in clinical settings.

Introduction

In the past few years, the topography of the cortical surface has raised a lot of interest, in particular to find biomarkers of pathologies [START_REF] Im | Quantification and discrimination of abnormal sulcal patterns in polymicrogyria[END_REF]; [START_REF] Auzias | Atypical sulcal anatomy in young children with autism spectrum disorder[END_REF] or to detect features associated with functional specificities (Sun et al. (2012b)). In order to perform automatic morphometry based on the local organization of the cortex, studying high-level objects such as cortical sulci allows to alleviate the dependency on the one-to-one voxel/vertex correspondence which is required in traditional Voxel-Based and Surface-Based Morphometry (VBM, [START_REF] Ashburner | Voxel-based morphometry-the methods[END_REF]; SBM [START_REF] Van Essen | An integrated software suite for surface-based analyses of cerebral cortex[END_REF]). This has led to the development of Object-Based Morphometry (OBM, [START_REF] Mangin | Object-based morphometry of the cerebral cortex[END_REF]), in which each sulcus is described by a large set of attributes. OBM has been successfully used to characterize various populations of subjects, for instance by mapping gender differences [START_REF] Duchesnay | Classification based on cortical folding patterns[END_REF] or by distinguishing patients from control subjects in pathologies such as schizophrenia [START_REF] Cachia | Cortical folding abnormalities in schizophrenia patients with resistant auditory hallucinations[END_REF] or autism [START_REF] Auzias | Atypical sulcal anatomy in young children with autism spectrum disorder[END_REF].

More recently, a specific attention has been brought to the deepest part of sulci, either to elaborate theoretical models of cortical anatomy and development [START_REF] Régis | Sulcal root generic model: a hypothesis to overcome the variability of the human cortex folding patterns[END_REF], or to automatically extract robust cortical landmarks [START_REF] Meng | Spatial distribution and longitudinal development of deep cortical sulcal landmarks in infants[END_REF]; [START_REF] Auzias | Deep sulcal landmarks: Algorithmic and conceptual improvements in the definition and extraction of sulcal pits[END_REF]. For the latter, the seminal work of [START_REF] Lohmann | Deep sulcal landmarks provide an organizing framework for human cortical folding[END_REF] has been particularly important in defining the concept of sulcal pits, with follow-ups brought by [START_REF] Im | Spatial distribution of deep sulcal landmarks and hemispherical asymmetry on the cortical surface[END_REF] and [START_REF] Auzias | Deep sulcal landmarks: Algorithmic and conceptual improvements in the definition and extraction of sulcal pits[END_REF]. This provides a finer spatial representation than the one offered by sulci -because a single sulcus can contain several sulcal pits. In particular, modelling local patterns of sulcal pits as graphs (see Fig. 2) has proved effective in order to establish links with genetic factors [START_REF] Lohmann | Deep sulcal landmarks provide an organizing framework for human cortical folding[END_REF]; [START_REF] Im | Quantitative comparison and analysis of sulcal patterns using sulcal graph matching: a twin study[END_REF] or characterize groups of patients [START_REF] Im | Quantification and discrimination of abnormal sulcal patterns in polymicrogyria[END_REF][START_REF] Im | Atypical sulcal pattern in children with developmental dyslexia and at-risk kindergarteners[END_REF]. These studies are all based on the statistical analysis of betweensubjects comparisons of pit-graphs defined within a pre-determined region of interest (ROI). These comparisons use a carefully designed similarity measure described in [START_REF] Im | Quantitative comparison and analysis of sulcal patterns using sulcal graph matching: a twin study[END_REF], that relies on a spectral matching method to define pit-to-pit correspondences between subjects.

The main limitations of this ROI-based method are the following. First, the effectiveness of the similarity measure of [START_REF] Im | Quantitative comparison and analysis of sulcal patterns using sulcal graph matching: a twin study[END_REF] strongly depends on the choice of the values of seven hyper-parameters which required a large set of specific experiments. Secondly, as with any ROI-based approach, one need to have strong a priori hypotheses to define the location and the size of the region to study. In that regard, the literature offers very little a priori information. The few published studies that looked at cortical folding patterns using sulcal pits examined either brain lobes in [START_REF] Im | Quantitative comparison and analysis of sulcal patterns using sulcal graph matching: a twin study[END_REF][START_REF] Im | Quantification and discrimination of abnormal sulcal patterns in polymicrogyria[END_REF] or single sulcal basins [START_REF] Im | Spatial distribution of deep sulcal landmarks and hemispherical asymmetry on the cortical surface[END_REF]; [START_REF] Auzias | Deep sulcal landmarks: Algorithmic and conceptual improvements in the definition and extraction of sulcal pits[END_REF]. Since both these approaches made it possible to detect significant effects at very different scales, it seems necessary to systematically conduct such study over a large range of spatial scales.

In order to overcome these limitations, we here introduce Structural Graph-Based Morphometry (SGBM), which is an extension of the work previously presented in [START_REF] Takerkart | Mapping cortical shape differences using a searchlight approach based on classification of sulcal pit graphs[END_REF]. This framework relies on three main contributions:

i) The design of a graph kernel that provides a new similarity measure between pit-graphs and allows -through the vast array of existing kernel methods [START_REF] Scholkopf | Learning with Kernels: Support Vector Machines, Regularization, Optimization, and Beyond[END_REF] -to perform various statistical analyses directly in graph space. This kernel has very few parameters that can be efficiently inferred from the data.

ii) The definition of a searchlight scheme -the first designed to perform brain morphometry -that yields information maps estimated from patterns of sulcal pits constructed at different spatial scales. Searchlight methods, introduced by [START_REF] Kriegeskorte | Information-based functional brain mapping[END_REF], consist in fitting a multivariate statistical model (e.g a classifier) on patterns defined in a local neighborhood, and repeating this operation in a sliding window fashion to fully cover the brain. A summary statistic (for instance the accuracy of the classifier) is then assigned to the center of each neighborhood, thus yielding a spatial information map that allows the localization of the informative regions.

iii) The construction of a non parametric multi-scale inference strategy that facilitates this localization by jointly analyzing the searchlight information maps obtained at all scales and offers a high detection power. Multi-scale methods aim at studying phenomena for which the optimal scale to be used is unknown [START_REF] Koenderink | The structure of images[END_REF]; [START_REF] Lindeberg | Scale-Space Theory in Computer Vision[END_REF], as when studying the local organization of sulcal pits. They have been used in neuroimaging for various tasks, such as the description of activation patterns in PET [START_REF] Coulon | Structural group analysis of functional activation maps[END_REF]) and fMRI [START_REF] Operto | Structural analysis of fMRI data: A surface-based framework for multi-subject studies[END_REF]) data or the segmentation of subcortical regions in anatomical MRI [START_REF] Wu | Hierarchical multi-atlas label fusion with multi-scale feature representation and label-specific patch partition[END_REF]). [START_REF] Etzel | Searchlight analysis: Promise, pitfalls, and potential[END_REF] also suggested that multi-scale strategies could be useful to desambiguate the a posteriori interpretation of regions detected by searchlight-based methods, which we implement here.

SGBM combines these three novelties in order to build a new fully automatic brain morphometry framework. This framework studies local patterns of sulcal pits by comparing them using our graph kernel, in order to define a classification-based searchlight scheme that overcomes the limitations of ROI-based approaches while embedding a multi-scale inference framework. This makes it possible to detect differences in local folding patterns at different scales between several groups of individuals. In the following, we first present the SGBM framework step by step, including a detailed description of our three contributions. We then demonstrate the power of our framework on two classical brain mapping problems for which complex patterns of anatomical differences have previously been reported in the literature: the mapping of asymmetries between the left and right brain hemispheres and the detection of cortical shape differences between male and female subjects.

Methods

Extracting sulcal pits from T1 MR images

In order to obtain the sulcal pits from an anatomical MR image, we first perform the segmentation of the cortical ribbon and the cortical reconstruction according to [START_REF] Dale | Cortical surface-based analysis: I. segmentation and surface reconstruction[END_REF] using the freesurfer image analysis suite1 . Then, we use the method describe in [START_REF] Auzias | Deep sulcal landmarks: Algorithmic and conceptual improvements in the definition and extraction of sulcal pits[END_REF] to extract the sulcal pits from the cortical sheet of each subject and the Depth Potential Function of [START_REF] Boucher | Depth potential function for folding pattern representation, registration and analysis[END_REF] to define a depth feature for each sulcal pit. This pit extraction method was carefully designed to obtain reproducible sulcal pits in every cortical region and not only in the deepest sulci, which is of critical importance to study patterns of pits centered around all cortical locations, as aimed at by the framework introduced in this paper. In order to compare pits between different subjects of the population S to be studied, we need to bring them into a common space. We use the cortical registration algorithm available in freesurfer [START_REF] Fischl | High-resolution intersubject averaging and a coordinate system for the cortical surface[END_REF] in order to align the convexity information across individuals after having inflated these meshes to a sphere. This sphere, onto which all subjects are aligned, will serve as our common space onto which we project the sulcal pits and basins of every subject, as illustrated on Fig. 1. For each subject s ∈ S, we thus obtain the set of sulcal pits P s = {π s i } i=N s i=1 localized on the sphere and their corresponding basins {β s i } i=N s i=1 . Note that the number of pits N s can vary across subjects.

Representing folding patterns as pit-graphs

A natural way to formally represent the local folding pattern of a given part of the cortex is to construct a region adjacency graph from the set of sulcal pits and basins located within this region, as proposed in [START_REF] Im | Quantitative comparison and analysis of sulcal patterns using sulcal graph matching: a twin study[END_REF]. Given a subset P ⊂ P s of N pits (N ≤ N s ), each pit π s i ∈ P defines a node of the graph. The graph edges are then given by the spatial adjacency of their associated basins: this defines a binary adjacency matrix A = (a ij ) ∈ R N ×N (a ij = 1 if β i and β j are adjacent, and 0 otherwise), that encodes the spatial organization of the folding pattern. We can then add attributes to each graph node in order to better characterize the folding pattern. In the present paper, we use two attributes for each pit: i) its euclidean coordinates on the common sphere X i , to compare graph nodes localization across subjects, and ii) its depth d i , because it is an intrinsic characteristic of the pit. The mapping from each subject to the common sphere explicitly minimizing metric distorsions [START_REF] Fischl | High-resolution intersubject averaging and a coordinate system for the cortical surface[END_REF], this sphere is a good space to compare pits locations across subjects in a spatially uniform manner. Note that other attributes can be added in the same manner without affecting the rest of the method, making these graphs a versatile representation. Let D = {d i } ∈ R N be the vector of depth values and X = {X i } ∈ R N ×3 be the matrix of coordinates of all graph nodes. Any pattern of sulcal pits can therefore be fully represented by an attributed graph defined as G = (P, A, D, X ). Examples of sulcal pit graphs are shown on Fig. 2.

Measuring pattern similarity using a graph kernel

Our first contribution is to introduce a graph kernel that provides a new similarity measure between pit-graphs. It is adapted from of a graph kernel that we designed to handle inter-subject variability in fMRI data [START_REF] Takerkart | Graph-based inter-subject pattern analysis of fMRI data[END_REF]. In order to compare two pit-graphs G = (P G , A G , D G , X G ) and H = (P H , A H , D H , X H ), our kernel exploits the different characteristics of the graphs by comparing all pairs of nodes (i.e potential edges) g ij and h kl , respectively in G and H. As such, it belongs to the class of walk-based graph kernels [START_REF] Gärtner | On graph kernels: Hardness results and efficient alternatives[END_REF] and uses the most elementary walks, of length one. It combines (see eq. 4) the different features of the graphs by using several sub-kernels within the convolution kernel framework [START_REF] Haussler | Convolution kernels on discrete structures[END_REF]. A first sub-kernel acts on the graph structure itself, by ensuring that the comparisons are performed only if g ij and h kl are actual edges. This is done with the linear kernel on the binary entries of the adjacency matrices:

k a (g ij , h kl ) = a G ij .a H kl , (1)
which takes the value 1 if a G ij = a H kl = 1 and 0 otherwise. A second sub-kernel uses a product of Gaussian kernels on the coordinates of the nodes of g ij and h kl to compare the locations of edges across subjects:

k x (g ij , h kl ) = e -X G i -X H k 2 /2σ 2 x • e -X G j -X H l 2 /2σ 2 x .
(2)

In practice, k x acts as a spatial filter that weights the comparisons of edges with their proximity, thus eliminating the comparisons of edges that are far away from eachother while tolerating inter-subject variability (if edges are close, but not perfectly matched across subjects). Finally, the last sub-kernel compares the depth attributes using the same principle:

k d (g ij , h kl ) = e -d G i -d H k 2 /2σ 2 d • e -d G j -d H l 2 /2σ 2 d (3)
Other optional attributes that could be added to further characterize the sulcal pits would be treated in the same manner.

The full kernel is defined as the combination of all sub-kernels applied on all pairs of nodes of G and H, which gives in our case:

K(G, H) = N G i,j=1 N H k,l=1 k a (g ij , h kl ) • k x (g ij , h kl ) • k d (g ij , h kl ) (4)
Note that the number of nodes N G and N H in G and H can be different. We then perform the following normalization

K(G, H) = K(G, H) K(G, G)K(H, H) , (5)
which ensures that ∀G, K(G, G) = 1, i.e that the similarity between a graph and itself does not depend on the graph. Note that this kernel has two hyper-parameters σ x and σ d , which are the bandwidths of the gaussian sub-kernels that respectively act on the coordinates and the depth features. In order to choose the values of these parameters, we use an extension of a standard practice whose effectiveness has been demonstrated in [START_REF] Takerkart | Graph-based inter-subject pattern analysis of fMRI data[END_REF]: we set σ x and σ d to the median euclidean distances between the coordinates and depth attributes of the nodes in all graphs available at training time. Note that in the case where the pit-graphs represent local folding patterns, as in the rest of the paper, the value of σ x will adapt to this property of the graphs and consequently, k x will act as a spatial filter within the neighborhood of interest. A: The sulcal pits and basins of a set of subjects projected onto the common spherical space. B: The Fibonacci point set used to sample the sphere. C: Around one of these points (larger green ball), the local pit-graphs defined for each subject within a spherical neiborhood (the node color encodes the depth of the sulcal pit). D: Classifying local pitgraphs into two categories (ex: patients vs. controls). E: Repeating this at all locations to obtain the searchlight information map: the Fibonacci points are projected onto the folded (left) or inflated (right) mesh (the sphere color encodes the amount of information available at this location for the classification task).

Searchlight mapping

Our second contribution is the definition of a searchlight method (Kriegeskorte et al. ( 2006)) that enables the mapping of cortical shape differences based on local folding patterns represented as sulcal pit-graphs. In order to construct a searchlight scheme for a given task, one needs to define the five following items: A: the domain of interest; B: the sampling strategy used to cover this domain; C: how to define local patterns at each location; D: the multivariate statistical model that addresses the task itself and E: the summary statistic to be used to create the information map. We now instantiate these five items, as illustrated on Fig. 3.

In order to define our pit-based searchlight scheme in a domain where the sulcal pits of all subjects can be compared, we use the common spherical space where all subjects have been aligned (see 2.1) as our domain of interest (item A). Then, the Fibonacci point set [START_REF] Niederreiter | Integration of nonperiodic functions of two variables by fibonacci lattice rules[END_REF]) provides a pseudo-evenly distributed sampling of the sphere, with an arbitrary number of points Q; this defines the sampling strategy of our searchlight scheme (item B). For a given location q ∈ {1, • • • , Q} and a given subject s ∈ S, we then define the set P s q,r of sulcal pits located within a radius r of point q, i.e P s q,r = {π s i ∈ P s | dist(q, π s i ) < r}, where dist designates the euclidean distance between two points. Given this set of pits, we can directly apply the graph construction scheme described in 2.2 to obtain the attributed graph that we note G s q,r . The set of graphs {G s q,r } s∈S defines the local patterns at location q (item C) that will serve as inputs to the statistical model defined hereafter.

Equipped with our graph kernel K defined in 2.3, we can use any kernel method to address a wide range of problems including regression, clustering or classification. Here, we focus on supervised classification to perform group discrimination and we use a non linear Support Vector Classifier as our multivariate statistical model (item D). Such a classifier is easy to estimate by using K, which is a positive definiete kernel, and the so-called kernel trickan idea that was introduced in [START_REF] Aizerman | Theoretical foundations of the potential function method in pattern recognition learning[END_REF]. We then assess the generalization power of the classifier by measuring its accuracy. In practice, for a given location q and a given value of r, we have a fully labeled dataset {(G s q,r , y s )} s=S s=1 at our disposal, where the label y s encodes the group subject s belongs to (for instance, y s = -1 for a control subject and y s = +1 for a patient). In order to obtain the accuracy ār q of the classifier, we resort to a 10-fold cross-validation which offers a good estimate of the true accuracy [START_REF] Kohavi | A study of cross-validation and bootstrap for accuracy estimation and model selection[END_REF].

The final step in setting up a searchlight scheme consists in defining a summary statistic from the output of the statistical model in order to construct the searchlight information map. Because the distribution of ār q is unknown, we first resort to a permutation scheme in order to estimate it under the null hypothesis of no differences between groups, as advocated in neuroimaging in general by [START_REF] Bullmore | Global, voxel, and cluster tests, by theory and permutation, for a difference between two groups of structural MR images of the brain[END_REF]; [START_REF] Nichols | Nonparametric permutation tests for functional neuroimaging: a primer with examples[END_REF] and more particularly in statistical analyses of searchlight information maps by [START_REF] Kriegeskorte | Information-based functional brain mapping[END_REF]; [START_REF] Stelzer | Statistical inference and multiple testing correction in classification-based multi-voxel pattern analysis (MVPA): Random permutations and cluster size control[END_REF]. Assuming the stationarity of this null distribution over the cortex, we pool together the classifier results at all searchlight locations to obtain an estimate of the cortex-wise distribution of ār q . This provides a map of uncorrected point-wise p-values

{p r q } q∈[1,••• ,Q]
. Then, at each point q, we compute the normal statistic that corresponds to p r q , which we note z r q ; this non linear transformation is in fact given by the so-called Inverse Error Function erf -1 , also called Normal Inverse Cumulative Distribution function. We use this z r q -score as our summary statistic (item E in the definition of our searchlight scheme), which defines the searchlight information map

Z r = {z r q } q∈[1,••• ,Q]
. This Z r map can easily be used to perform inference because we estimated its underlying distribution, and presents the nice property of being highly contrasted -as compared with the original accuracy map -thanks to the non linear transformation used to create the z r q -scores. The latter property will be useful at a later stage in our SGBM framework. The estimation of these information maps is implemented by the following algorithm, in which card designates the function that returns the size of the set given in input.

Algorithm 1: Computing the searchlight information map at scale r Input : a labeled dataset {(X s , y s )} s=S s=1

Input : a set of M permutations {Γ j } j=M j=1 , Γ 1 being the identity foreach permutation Γ j do foreach searchlight location q do compute ār q (j) from the re-labeled dataset {(G s q,r , y Γ j (s) )} s=S s=1 foreach permutation Γ j do foreach searchlight location q do p r q (j) = 1 Q×M card({(j ′ , q ′ ) | ār q ′ (j ′ ) ≥ ār q (j))} z r q (j) = erf -1 (p r q (j)) this defines the map Z r (j

) = {z r q (j)} q∈[1,••• ,Q] Output: a set of M maps {Z r (j)} j=M j=1
Output: the map corresponding to the true labels is Z r = Z r (1) Note that we also compute the set of information maps {Z r (j)} j=M j=1 obtained for all permutations because they will be used in the following. Also, by convention, the first permutation is set to be the identity, which means that the true information map Z r is included in this set as Z r (1).

Multi-scale spatial inference

In this section, we describe our third contribution, which is a multi-scale spatial inference strategy that both localizes the regions of the cortex where folding patterns are significantly different and estimates the spatial scale of these effects. First, we repeat the previously described searchlight procedure for a given set of ρ scales R . = {r 1 , • • • , r ρ }, yielding a set of information maps {Z r } r∈R , which each comes with its permuted versions {Z r (j)} j=M j=1 , forming a full set of M × ρ maps. Detecting at which locations q and scales r the classifier performance is above chance level -i.e testing z r q > 0, ∀(q, r) -is a massively univariate inference problem for which it is crucial to correct for multiple comparisons. Because these information maps exhibit correlations across both scales and space, an appealing way to deal with the multiple comparisons problem is to resort to cluster-based statistics, as introduced in [START_REF] Poline | Analysis of individual positron emission tomography activation maps by detection of high signal-to-noise-ratio pixel clusters[END_REF] and vastly used since in the neuroimaging literature. Following the non-parametric scheme described in [START_REF] Bullmore | Global, voxel, and cluster tests, by theory and permutation, for a difference between two groups of structural MR images of the brain[END_REF], this can be done in three steps:

❼ applying an arbitrary threshold to a statistical map to form contiguous regions of supra-threshold locations, classically referred to as clusters;

❼ defining a cluster-wise statistic and computing it for each cluster;

❼ estimating the null distribution of the maximum cluster statistic using a set of maps obtained by permuting the labels of the observations, which provides corrected cluster-wise p-values (see [START_REF] Nichols | Nonparametric permutation tests for functional neuroimaging: a primer with examples[END_REF]).

Thanks to the way our z r q statistic was designed (see 2.4), the initial cluster-forming operation is easy to perform by setting a threshold τ that corresponds to the desired one-sided uncorrected point-wise p-value; for instance, thresholding at z r q > τ with τ = 3.090 corresponds to p r q < 0.001. Note that because this is true for all values of r, a single threshold value can be used for all scales. Then, we need to select a cluster-wise statistic. The most commonly used one is the cluster size -the number of locations inside a cluster [START_REF] Poline | Analysis of individual positron emission tomography activation maps by detection of high signal-to-noise-ratio pixel clusters[END_REF], but it naturally favors large clusters. In order to avoid this bias, we use the cluster mass [START_REF] Bullmore | Global, voxel, and cluster tests, by theory and permutation, for a difference between two groups of structural MR images of the brain[END_REF]; [START_REF] Zhang | Cluster mass inference via random field theory[END_REF]), i.e the sum of supra-threshold point-wise statistics, which should better preserve the chances of small clusters to be detected thanks to the high contrast of our Z r maps. Note that this potential benefit brought by using the cluster mass -instead of the cluster size -would have been lesser if we had directly used the accuracy as our summary statistic. Finally, we introduce two complementary strategies to estimate the empirical distribution of the cluster-wise statistic in this multi-scale setting. They are described in the following and summarized on Fig. 4.

Single-scale clusters with corrections for multiple comparisons across

space and scales The first strategy consists in defining clusters by independently thresholding each Z r map. We therefore designate these as single-scale clusters in the rest of the paper. The multi-scale nature of this first inference method hails from the correction for multiple comparisons, which is performed in two steps. First, we estimate the null distributions of the maximum cluster mass separately for each single scale, by pooling the masses of clusters obtained with all permutations at the considered scale (see Alg. 2). This yields pvalues corrected for the multiple comparisons due to the repetition of tests in space. Note that it was necessary to distinguish null distributions between scales because they ended up being fairly different. In a second step, we used a Bonferroni correction to account for the repetition of tests across scales: multiplying the p-values obtained in the first step by the number of scales included in the analysis allows estimating corrected p-values that take into account the repetition of tests both in space and across scales, which can then be thresholded at the desired significance level to perform inference.

Algorithm 2: Computing spatially-corrected p-values for single-scale clusters Input : a scale r Input : a set of permuted single-scale maps {{Z r (j)} j=M j=1 , with the one obtained with the true labels for j = 1 Input : a threshold τ foreach permutation j do apply thresholding on Z r (j) with z r q (j) > τ get the set of clusters {c j i } i=N j i=1 from the thresholded map foreach cluster c j i do compute cluster mass m j i =

q∈c j i z r q (j) compute maximum mass M j = max i ({m j i } i=N j i=1 ) foreach true cluster c i = c 1 i do compute cluster mass m i = q∈c i z r q (1) compute corrected p-value p i = 1 M card({j | M j ≥ m i }) Output: single-scale clusters with their p-values {c i , p i } i=N 1 i=1
(corrected for multiple spatial testing) 2.5.2. Multi-scale clusters with corrections for multiple comparisons across space Our second strategy starts by designing a statistic that summarizes the results obtained across all scales at a given location q. Because we expect an effect to be observable across several consecutive scales, we first compute a sliding average of z r q across scales. Then, in order to detect effects potentially present at any scale, we select the maximum value from the results of the sliding average operation. This gives

Z R q = max r∈R ( 1 R r ′ , |r-r ′ |≤R/2 z r ′ q ), (6)
where R represents the size of the averaging window used across scales. This results in a single multi-scale information map

Z R = {Z R q } q∈[1,••• ,Q] ,
which summarizes the content of all single-scale maps {Z r } r∈R . We note {Z R (j)} j=M j=1 the set of corresponding permuted multi-scale maps. In the following, we will use the term multi-scale clusters to designate the clusters obtained by thresholding these multi-scale information maps. By changing the second input of Alg. 2 to this {Z R (j)} j=M j=1 set of permuted multi-scale maps, we can estimate the p-value -corrected for the repetition of tests performed in space -for each multi-scale cluster. Furthermore, at any location q within a significant cluster, we can estimate the preferred scale of the effect as the value of r that maximised the sliding average of z r q across scales for the true labels:

rR (q) = arg max r∈R ( 1 R r ′ , |r-r ′ |≤R/2 z r ′ q ) (7)
3. Experiments and results

Mapping gender and hemispheric differences

In order to explore the capabilities of our framework, we study two classical brain morphometry problems: the mapping of gender differences and cortical asymmetries. Because some neurological disorders are expressed differently in males and females (see [START_REF] Ruigrok | A meta-analysis of sex differences in human brain structure[END_REF]) and others are associated with abnormal cortical lateralization [START_REF] Toga | Mapping brain asymmetry[END_REF], studying these problems can directly provide information on the pathologies themselves. While they have already been studied using automatic morphometry methods such as VBM, SBM and OBM (see 3.3), our SGBM framework should shed a new light on both questions by examining differences in the organization of local folding patterns.

For this purpose, we use data from the Open Access Series of Imaging Studies2 (OASIS) cross-sectional database, which offers a collection of 416 subjects aged from 18 to 96. For each subject, three to four individual T1-weighted MP-RAGE scans were obtained on a 1.5T Vision system (Siemens, Erlangen, Germany) with the following protocol: in-plane resolution = 256x256 (1 mm x 1 mm), slice thickness = 1.25 mm, TR = 9.7 ms, TE = 4 ms, flip angle = 10➦, TI = 20 ms, TD = 200 ms. Images were coregistered and averaged to create a single image with a high contrast-to-noise ratio. From this database, we selected two groups of 67 male and 67 female healthy right-handed subjects, aged 18 to 34.

In practice, for the gender problem, we used the 134 subjects to set up a binary classification problem, with labels defined by y = -1 for a male subject and y = +1 for a female. For the asymmetry problem, an extra pre-processing step was carried out using the symmetric cortical template described in [START_REF] Greve | A surface-based analysis of language lateralization and cortical asymmetry[END_REF], as used in [START_REF] Auzias | Deep sulcal landmarks: Algorithmic and conceptual improvements in the definition and extraction of sulcal pits[END_REF] to study pit-based asymmetries. We then used the 134 symmetrized hemispheres of the male subjects to set up a binary classifcation problem where labels are y = -1 for a left hemisphere and y = +1 for a right hemisphere. For both problems, we use a Fibonacci point set of size Q = 2500, which provides a dense enough sampling of the sphere (see Fig. 3) to observe the expected continuity in the information maps estimated by the searchlight scheme (see results hereafter). The number of permutations was set to M = 5000, ensuring that the empirical distributions are accurately estimated. The crossvalidation used to estimate the average classification accuracy was carefully designed to ensure class balance in both the training and test datasets for each fold. As with any cluster-based analyses -which represent the vast majority of spatial inference methods used in neuroimaging, the initial clusterforming thresholding has an influence on the final results. With this in mind, we repeated our analysis with different threshold values and empirically selected the one that gave results which offered a good compromise between sensitivity and specificity. The results shown hereafter have been obtained with a cluster-forming threshold of τ = 3.090, which corresponds to a thresholding at p < 0.001 on the uncorrected point-wise p-values associated with the classification score. Finally, we used the following set of scales in order to study the multi-scale properties of local sulcal patterns: , 35, 40, 45, 50, 55, 60, 65, 70, 75, 80, 85, 90mm}. Each scale value defines the radius of a searchlight spherical neighborhood that is computed on the freesurfer sphere, i.e. the common space for all subjects, which has a 100mm radius. The smallest neighborhood (30mm Euclidean radius) is a spherical cap whose surface area is 2.25% of the sphere surface area. Equivalently, the largest neighborhood (90mm scale) represents 20.75% of the sphere area. In other words, we use a scale range that varies between 2.25% and 20.75% of an hemisphere surface area. This will yield patterns between the size of a small gyrus and one of a brain lobe, thus covering the smallest [START_REF] Im | Spatial distribution of deep sulcal landmarks and hemispherical asymmetry on the cortical surface[END_REF][START_REF] Auzias | Deep sulcal landmarks: Algorithmic and conceptual improvements in the definition and extraction of sulcal pits[END_REF] and largest [START_REF] Im | Quantitative comparison and analysis of sulcal patterns using sulcal graph matching: a twin study[END_REF][START_REF] Im | Quantification and discrimination of abnormal sulcal patterns in polymicrogyria[END_REF] effects previously reported in the literature when studying sulcal pits.

R . = {30

Results: methodological considerations

In this section, all the maps presented follow the same format. Each of the elements of the spherical Fibonacci point set is represented by a small sphere and projected onto the inflated cortical mesh. This allows for an easy visualization of the global maps.

Examining raw information maps

In order to assess the results provided by our SGBM framework, we first qualitatively observe raw single-scale and multi-scale information maps obtained with our pit-based searchlight scheme.

Single-scale maps Z r are shown on Fig. 5. They display some contrast at all scales r, i.e they contain locations of low and high z r q -values, sometimes reaching z r q = 5 or more -which corresponds to very significant uncorrected point-wise p-values (p r q < 10 -6 ). We can clearly observe the spatial correlation through the fact that the points of high statistic values are grouped together. Furthermore, the spatial smoothness of the maps and the size of the groups of searchlight locations with high z r q -values increase with r, which is a typical behavior of searchlight information maps [START_REF] Kriegeskorte | Information-based functional brain mapping[END_REF]. Finally, the persistence of the regions with high z r q -values across several consecutive values of r attests of the correlation across scales of the Z r maps, which also supports the relevance of the construction scheme used to define the multi-scale statistic Z R q (see Eq. 6). We now examine the multi-scale maps Z R , displayed on Fig. 6. Each of them summarizes the information contained at all scales r ∈ R. For R = 1, Z R q is the max of all z r q across all scales r; for R = card(R), i.e for R = 13 in our experiments, Z R q is the average of z r q across all scales. Overall, these maps are smoother than the single-scale Z r maps, which is expected because the construction of the multi-scale statistic Z R q acts as a smoothing operation across scales. For the smaller values of R, the contrast of the statistic map is weaker than the one offered by single-scale maps: some very large portions of the cortex show high values of Z R q , which means that the specificity offered by low values of R is limited. For very large values of R, the maps contain less locations with high statistical values, which is also expected because in this case, a high Z R q -value reflects a phenomena that exist over a larger number of consecutive scales. Therefore, it seems adequate to use an intermediary value of R, for which Z R q averages z r q over a limited number of consecutive scales.

Qualitative assessment of single-scale clusters

We then examine the single-scale clusters. Figs. 7 and8 show the clusters which present a p-value lower than 0.05, after correction for multiple comparisons across space and scales with the method described in Section 2.5.1. For both the gender and asymmetry problems, our framework detects significant clusters at all scales with r at least 40mm, demonstrating its detection power. Almost all clusters are persistent across scales, i.e for a given cluster at scale r, there exist clusters at nearby cortical locations for nearby scales. The clusters detected at finer scales are overall smaller than the ones detected with larger window sizes. This is caused by the larger smoothing effect of the large window sizes (as previously observed on the raw maps, see on Fig. 5); furthermore, we noted that the critical cluster mass corresponding to a corrected p-value of 0.05 (estimated by the procedure described in Section 2.5.1) increases with the size of the window, which therefore discards the smallest clusters at large scales. However, some small clusters remain significant at large scales. This overall demonstrates that our cluster-based spatial inference strategy is effective; in particular, the use of the cluster mass associated with our z r q statistic makes it possible to detect clusters at fine spatial scale and small clusters at larger scales, which is known to be challenging in practice in such a multi-scale context.

Qualitative assessment of multi-scale clusters

Figs. 9 and 10 show the significant clusters (corrected p < 0.05) obtained with the method described in 2.5.2, for R ∈ {1, 3, 5, 7, 9, 11, 13}. The observations made previously on the raw statistic map gets confirmed: some overly large clusters are detected for the smaller values of R, while nothing is detected for larger values (R ≥ 11). This confirms that an intermediary value of R might offer a satisfactory compromise between specificity and sensitivity. This is concordant with the intuition that an effect can live across several consecutive scales, but not over a very large number of consecutive scales as sought after with large R values.

Comparing single-scale and multi-scale results

We will now compare the clusters detected by our single-scale and multiscale approaches. Our objective is two-fold: i) to assess the consistency, or lack thereof, of the results of the two approaches, and ii) to facilitate the selection of the value of the R parameter of our multi-scale statistic.

For the gender problem, we examine Fig. 7 (single-scale clusters) and Fig. 9 (multi-scale clusters). The single-scale approach detects clusters in four main locations (tagged A, B, C and D on Fig. 7), which are present across a large number of scales. Their precise location and extent vary smoothly across scales. For instance, the cluster at location B slightly moves forward when the scale r increases; some merging phenomena also happen between clusters, as with the clusters at locations C and D, which are separated at scales r ≤ 60mm and fused into one larger cluster at scales r ≥ 70mm. The multi-scale approach detects clusters in similar locations for various values of the R parameter; however, the results with R = 7 are the only ones for which there are four separated clusters at the four locations A, B, C and D.

For the asymmetry problem, we examine Fig. 8 (single-scale clusters) and Fig. 10 (multi-scale clusters). The single-scale approach detects clusters in two main locations (tagged A and B on Fig. 8), which are present across almost the full range of scales from r = 30mm to r = 90mm. The clusters present around these locations also fuse together for r ≥ 80mm. The multiscale approach detects clusters in similar locations for various values of the R parameter; however, only the results obtained with R = 7 and R = 9 provide two separated clusters for each location A and B.

Overall, both our single-and multi-scale inference strategies offer satisfactory detection power. The fact that the significant clusters of both approaches occupy similar locations clearly demonstrates their consistency. In order to choose the value of the parameter R of the multi-scale strategy, we showed that choosing an intermediary value provides a good compromise between sensitivity and specificity, as found in our restults with R = 7 for the gender differences and asymmetry problems studied here.

Results: neuroscience considerations

In this section, we examine the neuroscientific relevance of the results obtained with our SGBM framework for the two problems of gender differences and cortical asymmetries, based on the multi-scale clusters detected with R = 7. Each map shown on Figs. 11 and 12 presents the set of searchlight locations (i.e the Fibonacci points displayed on the sphere on Fig. 3) included in a given significant cluster and projected onto the folded cortical mesh. This allows for a precise localization of the extent of each cluster on the cortex which can facilitate the interpretation.

Gender differences

Gender differences have mostly been studied using standard morphometry tools such as VBM, with notably a meta-analysis presented in [START_REF] Ruigrok | A meta-analysis of sex differences in human brain structure[END_REF] and SBM, for instance in [START_REF] Im | Gender difference analysis of cortical thickness in healthy young adults with surface-based methods[END_REF]; [START_REF] Lv | Gender consistency and difference in healthy adults revealed by cortical thickness[END_REF]). To our knowledge, the only sulcus-based study that used OBM is [START_REF] Duchesnay | Classification based on cortical folding patterns[END_REF].

Our SGBM framework detected four clusters that contain the center locations of searchlight windows in which the patterns of sulcal pits differ between genders, which are shown on Fig. 11. Cluster A is centered in the collateral sulcus of the left hemisphere. It is a fairly small cluster comprising 35 searchlight points, with an average preferred scale of r = 60mm (following the definition of r in 2.5.2). Its localization is consistent with regions that were significant in [START_REF] Ruigrok | A meta-analysis of sex differences in human brain structure[END_REF]; [START_REF] Im | Spatial distribution of deep sulcal landmarks and hemispherical asymmetry on the cortical surface[END_REF]; [START_REF] Lv | Gender consistency and difference in healthy adults revealed by cortical thickness[END_REF]. Cluster B is located in the right hemisphere, centered around the lateral occipital sulcus. It comprises 29 searchlight locations, has an average preferred scale of r = 65mm, and is concordant with very small spots appearing in the surface based studies only [START_REF] Im | Spatial distribution of deep sulcal landmarks and hemispherical asymmetry on the cortical surface[END_REF]; [START_REF] Lv | Gender consistency and difference in healthy adults revealed by cortical thickness[END_REF]). Moreover, the depth difference observed in the lingual sulcus in [START_REF] Duchesnay | Classification based on cortical folding patterns[END_REF] probably contributes to the effect detected in this cluster -which, in its single scale version at r = 80mm, includes the lingual sulcus on the internal face (see Fig. 7). Cluster C, located in the right hemisphere, has its center of mass in the cingulate gyrus. It is a large cluster with 98 searchlight locations with an average preferred scale of r = 75mm. Finally, cluster D is a large cluster centered in the middle frontal gyrus, comprising 130 searchlight locations and with an average preferred scale of r = 75mm. Small regions located within the large clusters C and D were also found significant in [START_REF] Ruigrok | A meta-analysis of sex differences in human brain structure[END_REF]; [START_REF] Im | Spatial distribution of deep sulcal landmarks and hemispherical asymmetry on the cortical surface[END_REF]; [START_REF] Lv | Gender consistency and difference in healthy adults revealed by cortical thickness[END_REF], demonstrating the overall consistency of our results with the literature.

Asymmetries

Cortical asymmetries have been studied using a large variety of methods, from VBM (see [START_REF] Good | Cerebral asymmetry and the effects of sex and handedness on brain structure: A voxel-based morphometric analysis of 465 normal adult human brains[END_REF] for instance) and SBM [START_REF] Van Essen | Parcellations and hemispheric asymmetries of human cerebral cortex analyzed on surface-based atlases[END_REF]) to OBM (see [START_REF] Duchesnay | Classification based on cortical folding patterns[END_REF]) and sulcal pits characterization (in [START_REF] Im | Spatial distribution of deep sulcal landmarks and hemispherical asymmetry on the cortical surface[END_REF]; [START_REF] Auzias | Deep sulcal landmarks: Algorithmic and conceptual improvements in the definition and extraction of sulcal pits[END_REF]).

By studying the organization of local patterns of sulcal pits, our SGBM framework was able to detect two clusters containing the center locations of searchlight windows in which cortical folding patterns present local asymmetries, which are shown on Fig. 12. The first one (Cluster A) is a large cluster comprising 295 searchlight locations, with an average preferred scale of r = 70mm, which is centered around the planum temporale and also encompasses the posterior part of the sylvian fissure, the supra-marginal gyrus and some parts of the superior temporal gyrus and sulcus. These regions are all individually reported in the literature dedicated to cortical asymmetries - [START_REF] Van Essen | Parcellations and hemispheric asymmetries of human cerebral cortex analyzed on surface-based atlases[END_REF]; [START_REF] Im | Spatial distribution of deep sulcal landmarks and hemispherical asymmetry on the cortical surface[END_REF]; [START_REF] Auzias | Deep sulcal landmarks: Algorithmic and conceptual improvements in the definition and extraction of sulcal pits[END_REF], all using univariate analyses, and [START_REF] Duchesnay | Classification based on cortical folding patterns[END_REF], which found a depth difference in the superior temporal sulcus. This shows that the multivariate nature of our analysis allows to detect such locally distributed effects among a complex pattern of sulcal basins. Cluster B comprises 205 searchlight locations centered around the frontal pole and presents an average preferred scale of r = 60mm. It includes part of the lateral prefrontal cortex, and extends on the internal face to the middle part of the cingulate sulcus. These regions have been consistently detected with morphometry approaches such as VBM [START_REF] Good | Cerebral asymmetry and the effects of sex and handedness on brain structure: A voxel-based morphometric analysis of 465 normal adult human brains[END_REF]) and also with the study of sulcal pits frequency described in [START_REF] Auzias | Deep sulcal landmarks: Algorithmic and conceptual improvements in the definition and extraction of sulcal pits[END_REF] -which used the same pit extraction algorithm as ours. These concordant results corresponds to the frontal petalia, a protrusion of the anterior part of left hemisphere towards the right one [START_REF] Toga | Mapping brain asymmetry[END_REF].

Discussion

We have introduced Structural Graph-Based Morphometry (SGBM), a new fully automatic brain morphometry framework based on a multi-scale searchlight scheme that enables to localize differences in graphs of sulcal pits between populations. This is the first time that a searchlight scheme is designed for anatomical MR images of the brain and embedded in a multi-scale inference strategy. The main strengths of our framework are i) to capitalize on structural pattern recognition tools in order to discriminate complex cortical folding patterns represented as graphs of sulcal pits, ii) to overcome the intrinsic limitations of approaches that focus on a pre-defined region thanks to the searchlight approach, iii) to embed this searchlight approach into a powerful cluster-based multi-scale inference strategy.

The results obtained in this paper demonstrate the detection power of our SGBM framework in order to localize regions where sulcal pit-graphs are significantly different between groups, thus showing its ability to overcome the intrinsic limitations of previous ROI-based approaches. The multi-scale character of our SGBM framework was called for by the lack of a priori knowledge on the size of the relevant patterns. The multi-scale nature of these effects is clearly confirmed by our results: on the one hand, single-scale clusters are detected over a large span of scales, with varying locations and sizes (see Figs. 7 and8), and on the other hand, the detection of multi-scale clusters (see Figs. 9 and 10) shows that these effects are present over several consecutive scales. These multi-scale clusters also come with an estimate of the preferred scale of the effect, which vary between 60 and 75mm for the studied problems. These values, which are defined in freesurfer's common spherical space, correspond to intermediary region sizes between the smallest and largest previously reported effects, i.e to folding patterns comprising several sulci. This confirms that local interactions between neighboring sulcal basins or sulci can be informative. Moreover, our two inference strategies showed their complementarity. The multi-scale scheme is computationally efficient and offers a compact representation of the results, summarized in a single map, which eases the interpretation of the results. Comparatively, the single-scale clusters provide a more detailed scale by scale description of the phenomena at hand, but come at a higher computational cost -because the number of permutations to perform is multiplied by the number of scales considered. And because of the large number of resulting maps -one per scale, the results of the single-scale scheme might be more difficult to in-terpret. Future work could address this latter drawback by modelling more precisely the relations between clusters obtained at different scales, as was done in [START_REF] Operto | Structural analysis of fMRI data: A surface-based framework for multi-subject studies[END_REF] for fMRI data. Finally, note that our multi-scale inference strategy could benefit to other searchlight applications, in particular for fMRI data analysis where it has grew very popular; indeed, [START_REF] Etzel | Searchlight analysis: Promise, pitfalls, and potential[END_REF] suggested to use a systematic multi-scale approach to provide more robust statistical results and thus improve the interpretability, and our method directly fullfills this need.

The fact that the results obtained with our SGBM framework are in line with the existing literature also validate the non linear Support Vector Classifier as an adequate statistical model to study patterns of sulcal pits and their differences between populations. This demonstrates that the graph kernel that we introduced in 2.3 is a valid similarity measure for comparing pit-graphs, thus providing an alternative to the one of [START_REF] Im | Quantitative comparison and analysis of sulcal patterns using sulcal graph matching: a twin study[END_REF]. It is to be noted that the choice of the weights that define Im's metric required a set of dedicated experiments (see [START_REF] Im | Quantitative comparison and analysis of sulcal patterns using sulcal graph matching: a twin study[END_REF]) and the result of their followup work [START_REF] Im | Quantification and discrimination of abnormal sulcal patterns in polymicrogyria[END_REF][START_REF] Im | Atypical sulcal pattern in children with developmental dyslexia and at-risk kindergarteners[END_REF]) heavily depends on the values of these seven hyper-parameters. In contrast, our kernel only has two hyper-parameters that we directly estimate from the data. The only available element of comparison between the two measures is their ability to detect gender differences. In [START_REF] Im | Quantification and discrimination of abnormal sulcal patterns in polymicrogyria[END_REF], the authors reported no gender differences in their control group (first sentence in Results, p.3011), whereas our technique was able to find several regions with highly significant differences between male and female subjects. These differences can be explained by several reasons: i) our dataset was much larger (134 vs. 26 subjects); ii) our pit detection algorithm is slightly different; iii) their study focused on predetermined large regions (brain lobes), whereas our framework automatically estimates the localization and the size of the region; iv) our inference strategy is more powerful, thanks to the exploitation of the local redundancy that exist between locations within a cluster. In theory, it is also possible that our kernel provides a more expressive graph similarity measure, but we have no element to concur with this possibility. On the contrary, it was shown -respectively in [START_REF] Im | Spatial distribution of deep sulcal landmarks and hemispherical asymmetry on the cortical surface[END_REF] and in [START_REF] Takerkart | Graph-based inter-subject pattern analysis of fMRI data[END_REF] (where a kernel with a similar design was used on fMRI patterns) -that both similarity measures effectively use all elements of the graphs -their structure, the location of the graph nodes and their attributes. We therefore believe that they are sensitive to the same features of the folding pattern, that include: ❼ differences in sulcal depth distributed across a local folding pattern, as suggested by the concordant results with [START_REF] Duchesnay | Classification based on cortical folding patterns[END_REF] which detected depth differences co-localized with some of the SGBM clusters,

❼ differences in the local number of pits, as showed by the consistent results of SGBM and the pit frequency study of [START_REF] Auzias | Deep sulcal landmarks: Algorithmic and conceptual improvements in the definition and extraction of sulcal pits[END_REF],

❼ differences in pits localization, as suggested by the agreement between our results and the one of [START_REF] Im | Spatial distribution of deep sulcal landmarks and hemispherical asymmetry on the cortical surface[END_REF] which examined pits displacements;

❼ local organization differences, which might be induced by the two previous types of effects.

When it comes to interpreting the results provided by our SGBM framework, we demonstrated that a visualization tool such as the one used in 3.3 allows to easily localize the detected clusters on the cortex, which makes it possible to directly compare SGBM with other automatic morphometry methods such as VBM, SBM and OBM. For the gender differences and asymmetry problems, our results are consistent with previously reported differences in grey matter volume (VBM), cortical thickness (SBM) and sulcal shape (OBM), which was expected because these different features are -at a certain level -linked together. But for the first time, SGBM makes it possible to automatically detect markers associated with the spatial organization of cortical folding. To go beyond a simple reporting of the localization of the detected clusters and in order to unveil the neuroscientific underpinnings of such effect, there is a crucial need to understand which of the four aforementioned characteristics of the folding patterns actually contributes to the group discrimination in a given cluster. This is a very challenging problem because of the complex nature of cortical folding and because our graph kernel combines several feature types in a single similarity measure. But its modular design makes it possible to conduct posthoc sensitivity analyses similarly to what was done in [START_REF] Im | Quantitative comparison and analysis of sulcal patterns using sulcal graph matching: a twin study[END_REF][START_REF] Im | Quantification and discrimination of abnormal sulcal patterns in polymicrogyria[END_REF][START_REF] Im | Atypical sulcal pattern in children with developmental dyslexia and at-risk kindergarteners[END_REF] and in [START_REF] Takerkart | Graph-based inter-subject pattern analysis of fMRI data[END_REF], which simply consists in adding or removing a feature, repeating the analysis and assessing its contribution by examining the performance difference. Finally, because our similarity measure is a positive definite kernel, we can exploit the vast array of available kernel methods to perform other types of posthoc analyses that could help us gain further neuroscientific insights.

For instance, we could perform dimensionality reduction using kernel PCA in order to examine the dimensions in graph-space that explain the largest variance in a given cluster. Then, similarly to what was done in Sun et al. (2012a) using another dimensionality reduction technique, it is possible to project the pit-graphs of all subjects along each principal dimension and to plot them along the corresponding axis in order to visually assess the main factor(s) of influence. Beyond these two leads -sensitivity analyses and dimentionality reduction, the intrinsic complexity of the folding pattern leads us to believe that further work will be needed to better understand what contributes to the differences detected with our SGBM framework and how.

Conclusion

We have introduced a new brain mapping technique dedicated to studying differences in the local organization of cortical anatomy, by the mean of patterns constructed from the deepest part of the cortical sulci -the sulcal pits. Our technique -that we called Structural Graph-Based Morphometry (SGBM), is the first morphometry framework that combines a graph kernel to compare sulcal pit-graphs, a searchlight scheme to overcome the limitations of ROI-based approaches and a multi-scale inference strategy that allows to find significant effects at different spatial scales. We have demonstrated that SGBM offers good detection power on two problems that respectively examined the anatomical gender differences and the cortical asymmetries. This versatile and powerful pit-based morphometry framework should therefore find numerous applications, in particular to study cortical development and clinical populations, two fields where examining the organization of sulcal patterns at different spatial scales is particularly relevant.

Vitae

Sylvain Takerkart is a research engineer at the French Centre National de la Recherche Scientifique (CNRS), working as the manager of the scientific computing facility of the Institut de Neurosciences de la Timone. His research activities focus on developing and applying machine learning tools for new image processing problems encountered in neuroscience.

Guillaume Auzias received the Masters degree in Applied Mathematics from University Pierre and Marie Curie -Paris 6 in 2006 and his PhD at the Brain and Spine (ICM) institute in 2009. The main goal of his work is to demonstrate the benefit of integrating information derived from cortical folds in order to define correspondences between subjects.

Lucile Brun is a research engineer at the Institut de Neurosciences de la Timone. With a dual qualification in medical imaging and neuroscience, she is a member of the Social Cognition across Lifespan and Pathologies research group, which studies morphology changes in neuro-developmental pathologies. She contributes to the design of extraction methods for new anatomical features such as the sulcal pits.

Olivier Coulon is a director of research at the Laboratoire des Sciences de l'Information et des Systèmes, CNRS, Aix-Marseille University, and the head of the Methods and Computational Anatomy research group (http: //www.meca-brain.org). For more than 15 years his research has been dedicated to neuroimaging data analysis. He is particularly interested in modelling and quantifying cortical organization and variability, and the link between anatomy and function.

Figure 1 :

 1 Figure 1: Left: sulcal pits of one subject, represented on the cortical sheet. Right: after cortical registration to the template and projection to the common spherical space, in which inter-subject pits locations are compared.

Figure 2 :

 2 Figure 2: Sulcal pit-graphs in the right frontal lobe for four subjects, displayed on the inflated cortex. The node color encodes the depth of the pit, from blue (shallow) to red (deep).

Figure 3 :

 3 Figure 3: Illustrations of the pit-based searchlight scheme used in our SGBM framework.A: The sulcal pits and basins of a set of subjects projected onto the common spherical space. B: The Fibonacci point set used to sample the sphere. C: Around one of these points (larger green ball), the local pit-graphs defined for each subject within a spherical neiborhood (the node color encodes the depth of the sulcal pit). D: Classifying local pitgraphs into two categories (ex: patients vs. controls). E: Repeating this at all locations to obtain the searchlight information map: the Fibonacci points are projected onto the folded (left) or inflated (right) mesh (the sphere color encodes the amount of information available at this location for the classification task).

Figure 4 :

 4 Figure 4: Illustrations of the two inference strategies. A: the set of permutations of the labels (each column is a permutation, black corresponds to y = -1 and white to y = +1; the first permutation in the left-most column contains the true labels). B: single-scale searchlight information maps, for each scale and each permutation. C & D: single scale clusters, and the null distributions of their maximum cluster mass, estimated from the set of permuted single-scale clusters at each scale. E, F & G: multi-scale information maps, multi-scale clusters, and the null distribution of their maximum cluster mass. The distribution shown in D and G are then used to assess the corrected p-values, respectively for the single-scale and the multi-scale clusters obtained with the true labels.

Figure 5 :

 5 Figure 5: Examples of raw single-scale information maps Z r for different neighborhood radiuses r (gender differences, external face of the right hemisphere).

Figure 6 :

 6 Figure 6: Examples of raw multi-scale information maps Z R (gender differences, external face of the right hemisphere). Each of these maps summarizes the information from all single-scale Z r maps shown on Fig 5. This multi-scale statistic depends on the number R of consecutive scales r over which the single-scale statistic is averaged (see Eq. 6).

Figure 7 :

 7 Figure 7: Gender differences. For each value of r, four views are shown, two per hemisphere: the left hemisphere is shown on the left, and the right on the right, with respectively their external and internal faces at the top and the bottom. The significant single-scale clusters (corrected p < 0.05) are colored with their p-value. Four main clusters, tagged A, B, C and D, appear to be persistent across scales, sometimes after having merged.

Figure 8 :

 8 Figure 8: Asymmetries. For each value of r, the single-scale clusters presenting significant differences between the left and right hemispheres (corrected p < 0.05) are shown on the external and internal faces of the left hemisphere. Two main clusters, tagged A and B, appear to be persistent across scales.

Figure 9 :

 9 Figure 9: Gender differences. Simlarly to Fig.7, four views are presented for each value of R, showing thresholded significant multi-scale clusters (corrected p < 0.05). For R = 7, four clusters survive, at similar locations as single-scale ones; they are tagged A, B, C and D, accordingly to Fig. 7.

Figure 10 :

 10 Figure 10: Asymmetries. Significant muli-scale thresholded clusters. For R = 7, the same two clusters as in the single-scale analysis are detected; they are tagged A and B, accordingly to Fig. 8.

Figure 11 :

 11 Figure 11: Gender clusters. Each map presents the Fibonacci points included in a given cluster, among the four significant clusters detected using our multi-scale statistic Z R and R = 7. The larger green sphere represents the center of mass of the cluster.

Figure 12 :

 12 Figure 12: Asymmetry clusters. These maps each present one of the two significant clusters detected using our multi-scale statistic Z R with R = 7. The larger green sphere represents the center of mass of the cluster. For clarity, the map of Cluster B (right) is presented under two different views.

www.freesurfer-software.org

www.oasis-brains.org

Acknowledgments

We would like to thank the anonymous reviewers for their constructive suggestions and comments. G. Auzias and L. Brun were funded through research grants from Fondation de France (OTP 38872) and Fondation Orange (S1 2013-050). The research of S. Takerkart and O. Coulon is supported by the French Centre National de la Recherche Scientifique.