Pietro Gori

Olivier Colliot

Linda Marrakchi-Kacem

Yulia Worbe

Fabrizio De

Vico Fallani

Mario Chavez

Cyril Poupon

Linda Marrakchi

A Bayesian Framework for Joint Morphometry of Surface and Curve meshes in Multi-Object Complexes. submitted to Medical Image Analysis. Chapter 2 •

to be associated to a dierent wiring architecture of the white matter with respect to healthy subjects. This would involve dysfunctions of the cortico-striatothalamo-cortical (CSTC) circuits of the brain which could underlie the pathological symptoms.

The CSTC neural loops, schematically represented in Fig. 1.1, are topographically and functionally well separated into sensorimotor, associative and limbic circuits which are implicated in motor, cognitive and motivational aspect of behaviour respectively [START_REF] Alexander | Parallel Organization of Functionally Segregated Circuits Linking Basal Ganglia and Cortex[END_REF]. Every circuit is composed of neural projections from several functionally related cortical areas which are sent to restricted portion Chapter 1. Introduction of the striatum (caudate nucleus and putamen). These striatal regions send further projections to the globus pallidus and substantia nigra, which in turn project to a specic region of the thalamus. Each thalamic region projects back to one (or more) cortical areas. The role of the basal ganglia and thalamus is to integrate the information coming from dierent inputs of the entire cerebral cortex and to "funnel" back these inuences again to particular areas of the cortex. Dysfunctions in these circuits are thought to be associated with organisational and morphological alterations of their anatomical components. The analysis of the conguration of these structures, namely their relative position, implies a systemic approach. All the components of the circuits, from both white and grey matter, need to be studied as a whole and not as a set of independent components. Moreover, neurodevelopmental disorders are characterised by an important clinical heterogeneity. Each symptom is thought to be related to a specic pattern of anatomical alterations [START_REF] Tuite | Magnetic Resonance Imaging in Movement Disorders: A Guide for Clinicians and Scientists[END_REF]). Thus, it seems more likely that a syndrome would present a range of abnormalities rather than a single specic one. This pathological variability would be added to the normal variability presents in the healthy population.

The analysis of neurodevelopmental disorders need therefore holistic and statistical approaches to estimate the typical brain conguration associated to a disorder and its variability. This should then be compared with the characteristic organisation of a population of healthy controls in order to distinguish the normal variations 1.2. Gilles de la Tourette syndrome from the pathological ones. Alterations should be clearly localised, quantied and characterised in order to be easily interpreted by clinicians and permit to better characterise the pathophysiological mechanisms underlying the disorder. Furthermore, the estimated typical organisations of both patients and healthy controls, together with their variability, could be employed to classify the conguration of a new patient in order to improve the diagnostic accuracy.

To this end, we propose in this thesis a computational and statistical framework to analyse morphological and organisational anomalies altering the anatomy of the neural circuits in neurodevelopmental disorders and we evaluate it on a disease model about Gilles de la Tourette syndrome.

Gilles de la Tourette syndrome

Gilles de la Tourette syndrome (GTS) is a heritable, childhood onset, neuropsychiatric disorder characterized by the presence of sudden, rapid, recurrent and nonrhythmic movements (motor tics) or vocalisations (phonic tics) [START_REF] Hallett | Tourette Syndrome: Update[END_REF]. In more than 80-90% of patients, tics are accompanied by psychiatric comorbidities such as obsessive-compulsive disorder (OCD) [START_REF] Kimber | An Update on Tourette Syndrome[END_REF]. Typically the disorders begin in early childhood (6-7 years) with transient bouts of simple motor tics which are usually followed several years later by phonic tics. Tics may initially wax and wane but eventually they become persistent reaching the maximum severity at an age around 10-12 years. After that there is usually a steady decline in symptoms [START_REF] Leckman | Tourette's syndrome[END_REF][START_REF] Jankovic | Tourette syndrome: Evolving concepts[END_REF]. Once thought to be rare, current estimates change this view. GTS occurs in all countries, ethnic groups and social classes and epidemiological studies suggest that its overall prevalence could be as high as 1% of the general paediatric population in most countries [START_REF] Leckman | Tourette's syndrome[END_REF][START_REF] Stern | Gilles de la Tourette's syndrome and its impact in the UK[END_REF][START_REF] Jankovic | Tourette syndrome: Evolving concepts[END_REF].

The exact aetiology of GTS is still unclear as well as its underlying pathophysiological mechanisms [START_REF] Worbe | Neuroimaging of Gilles de la Tourette syndrome[END_REF]. Dierent living environment and lifestyle factors are thought to be implicated in the development of GTS (Mc-Naught and Mink, 2011). Genetic factors have also a major role and the clinical heterogeneity of this syndrome suggests that it might be related to dierent genes and to several genetic variants. Unfortunately, at present the genes that predispose individuals to develop the disorder are unknown [START_REF] Leckman | Tourette's syndrome[END_REF].

Driving hypothesis about GTS pathophysiology

In the last years, there has been an increasing interest for the neuroanatomical causes of this syndrome due also to the arrival of new high-resolution imaging tools.

Dierent hypotheses have been proposed. Among them, inference from various approaches support that GTS is associated with dysfunctions of the cortico-striatothalamo-cortical (CSTC) circuits [START_REF] Worbe | Altered structural connectivity of cortico-striatopallido-thalamic networks in Gilles de la Tourette syndrome[END_REF][START_REF] Singer | Motor control, habits, complex motor stereotypies, and Tourette syndrome[END_REF][START_REF] Ganos | The functional anatomy of Gilles de la Tourette syndrome[END_REF][START_REF] Mcnaught | Advances in understanding and treatment of Tourette syndrome[END_REF]. According to that, GTS symptoms would arise Chapter 1. Introduction due to an atypical integration of information which would be associated to anatomical and functional abnormalities in the CSTC circuits. Alterations could be the result of a dierent brain maturation in GTS patients. In a typical brain development there is an initial increase in the number of neural connections which is followed by a selective pruning during brain maturation which optimizes the eciency of the neural connections. An erroneous pruning in GTS, together with defects in the axonal myelinization, could lead to less ecient CSTC circuits and consequently also to the development of the symptoms. The heterogeneous clinical expression of GTS could be explained by dysfunctions impairing dierent CSTC circuits. For instance, simple tics are thought to be associated with abnormalities in the motor circuit [START_REF] Worbe | Distinct structural changes underpin clinical phenotypes in patients with Gilles de la Tourette syndrome[END_REF][START_REF] Fahim | Somatosensorymotor bodily representation cortical thinning in Tourette: Eects of tic severity, age and gender[END_REF][START_REF] Sowell | Thinning of sensorimotor cortices in children with Tourette syndrome[END_REF] whereas complex tics could be linked to dysfunctions of the associative circuit [START_REF] Worbe | Distinct structural changes underpin clinical phenotypes in patients with Gilles de la Tourette syndrome[END_REF][START_REF] Neuner | White-matter abnormalities in Tourette syndrome extend beyond motor pathways[END_REF] and the presence of psychiatric co-morbidities, such as OCD and ADHD, could instead be associated with abnormalities in the limbic circuit [START_REF] Peterson | Morphologic Features of the Amygdala and Hippocampus in Children and Adults With Tourette Syndrome[END_REF][START_REF] Worbe | Distinct structural changes underpin clinical phenotypes in patients with Gilles de la Tourette syndrome[END_REF]. These hypotheses were corroborated by experiments based on primate models [START_REF] Worbe | Behavioral and Movement Disorders Induced by Local Inhibitory Dysfunction in Primate Striatum[END_REF](Worbe et al., , 2013)). Micro-injections of bicuculline, a γ-aminobutyric acid (GABA)ergic antagonist, were performed in dierent locations of the CSTC circuits to a group of 8 monkeys. These injections induced neuronal perturbations similar to three GTS symptoms, namely tic-like movements, hyperactivity and stereotyped behaviours. To localize the neurons perturbed by the bicuculline, the authors performed injections of axonal tracers followed by an histological analysis. It was shown that animals with abnormal movements resembling to simple motor tics had dysfunctions in the premotor and sensorimotor circuits. Instead, monkeys presenting behavioural disorders resembling to complex tics and compulsions had dysfunctions in the associative and limbic circuits.

These results were consistent with the aforementioned hypothesis but, as pointed out by the authors, they can not be directly translated to humans. Post-mortem studies in humans, like in [START_REF] Kataoka | Decreased number of parvalbumin and cholinergic interneurons in the striatum of individuals with Tourette syndrome[END_REF], have also corroborated this hypothesis. Unfortunately, the small sample size of brain specimens, the large clinical heterogeneity of GTS, the presence of comorbidities and its evolution with age can confound the clinical conclusions. Nowadays, the precise mechanisms of the dysfunctions acting on the CSTC circuits are still poorly understood.

For this reason, the role of neuroimaging has already become fundamental since it is a non-invasive and painless technique which could help answering fundamental questions such as Which are the anatomical abnormalities associated to the dysfunction of the CSTC circuits? and more importantly Where are they located and which is their extent? 1.4 The role of structural and diusion imaging Structural and diusion imaging studies seek to characterise the anatomy of the brain and the abnormalities associated to a pathology.

The role of structural and diusion imaging

Magnetic Resonance Imaging (MRI) is a radiologist technique used to visualize internal structures of the body in detail. It makes use of the property of nuclear magnetic resonance (NMR) to produce high quality two-or three-dimensional images. Dierent modalities exist and two of the most employed are the T1-weighted (T1-w) and the diusion MRI (dMRI). The rst one provides an excellent soft tissue contrast which is especially useful for characterising the anatomy of both grey and white matter in the brain. The second one provides instead images whose contrast is inuenced by the dierences in water molecule diusivity and it is used to study the micro-structure and the organisation of the white matter axons.

Dierent studies based on these two modalities have already produced a reasonable body of evidence consistent with the aforementioned hypothesis about GTS.

For instance, morphometric studies on T1-w scans have shown a volume reduction in the caudate and in the sensorimotor, parietal, orbitofrontal and anterior cingulate cortical areas. Moreover, it was also reported an enlarged hippocampus, amygdala, putamen and corpus callosum in GTS patients [START_REF] Felling | Neurobiology of Tourette Syndrome: Current Status and Need for Further Investigation[END_REF]. These results were obtained by comparing scans of healthy subjects (controls) with GTS patients using mainly two techniques. The rst one is based on the comparison of morphometric measurements, such as the volume, of single and clearly dened grey matter structures of the circuits (i.e. caudate and putamen) which are manually or automatically segmented [START_REF] Felling | Neurobiology of Tourette Syndrome: Current Status and Need for Further Investigation[END_REF][START_REF] Peterson | Basal ganglia volumes in patients with gilles de la tourette syndrome[END_REF]. This technique produces results which are rather easy to interpret but they give just an overall insight on the pathological anatomical eects of GTS without being able to reveal the exact location of the abnormalities and whether inter-linked neighbour components are also aected.

The second technique is called Voxel-Based Morphometry (VBM) [START_REF] Ashburner | Voxel-Based MorphometryThe Methods[END_REF] and it is based on non-rigid registrations of the structural scans of both controls and patients to a reference scan called template. This permits to bring all the scans in the same reference frame allowing a voxel-wise comparison between the deformed images. After segmenting every subject scan into dierent tissue types and smoothing the resulting segmentations, mass-univariate statistics are used to detect statistically signicant dierences in the local concentration of grey matter between the two groups. A similar approach, called Tensor-Based Morphometry (TBM) [START_REF] Davatzikos | A computerized approach for morphological analysis of the corpus callosum[END_REF], uses instead the deformations to quantify and model the anatomical dierences between subjects. It is based on the comparison of the Jacobian matrix elds of the deformations (or their determinants) using multi-variate (or mass univariate) statistics. These two techniques can also be combined in an augmented version of VBM where every voxel of the deformed image is multiplied by its volume change which is encoded in the Jacobian determinant of the deformation eld [START_REF] Ashburner | Voxel-Based MorphometryThe Methods[END_REF]. In this way it is possible to compare the absolute amount of gray matter instead than the relative concentration of grey matter with respect to the other tissues like white matter.

These methods do not focus on single structures and they result in a comprehensive evaluation of anatomical alterations across the whole brain. However, they can investigate only focal dierences at each voxel without taking into consideration Chapter 1. Introduction the organisation of the dierent components. Moreover, these techniques are highly sensitive to misalignement and misclassication errors which can bring to erroneous interpretations. Their target is usually the grey matter but they could be used to examine also the volume variations in the white matter. However, white matter in T1-w scans is usually characterized by large homogeneous regions with only subtle changes in intensity.

For this reason, white matter is predominantly analysed with dMRI scans where the diusion at each voxel is modelled as a tensor (Diusion Tensor Imaging -DTI) or using more accurate model-free approaches such as Q-Ball Imaging (QBI) [START_REF] Tuch | Q-ball imaging[END_REF]) which can resolve intra-voxel ber crossing, contrarily to DTI. White matter microstructure is then characterised by scalar quantities such as fractional anisotropy (FA) [START_REF] Pierpaoli | Toward a quantitative assessment of diusion anisotropy[END_REF] or generalized fractional anisotropy (GFA) [START_REF] Tuch | Q-ball imaging[END_REF][START_REF] Cohen-Adad | Detection of multiple pathways in the spinal cord using q-ball imaging[END_REF] which are based on DTI and QBI respectively. They dene the anisotropy level of diusion at each voxel and they can be interpreted as measures of the structural integrity of the ber bundles. These quantities can be compared voxel-wise across groups of subjects using VBM or TBSS (Tract-based spatial statistics) [START_REF] Smith | Tract-based spatial statistics: Voxelwise analysis of multi-subject diusion data[END_REF]. The latter method, similarly to VBM, starts by aligning all FA (or GFA) images to a template using a cascade of ane and non-linear registration. These images are then projected to a skeletonised mean FA (or GFA) image, which represents the centres of all tracts common to the group of subjects, where mass-univariate statistics are carried out.

Dierently from VBM, this technique permits to avoid the image smoothing step, where it is still unclear how to x the optimal smoothing parameter, and it is supposed to reduce the registration errors thanks to the projection to the mean FA (or GFA) image. Both techniques have been employed for studying GTS and with VBM it was shown a decrease in FA in the corticospinal tract, internal capsule and corpus callosum as well as an increase in the white matter underlying the ventral and posteriolateral part of the thalamus in patients subject to GTS [START_REF] Thomalla | Structural changes in the somatosensory system correlate with tic severity in Gilles de la Tourette syndrome[END_REF][START_REF] Neuner | White-matter abnormalities in Tourette syndrome extend beyond motor pathways[END_REF] whereas it was shown no dierences in FA using TBSS [START_REF] Govindan | Abnormal water diusivity in corticostriatal projections in children with Tourette syndrome[END_REF]. However, in another study based on TBSS [START_REF] Liu | Structural Abnormalities in Early Tourette Syndrome Children: A Combined Voxel-Based Morphometry and Tract-Based Spatial Statistics Study[END_REF], the authors found a signicant axial diusivity (AD) and mean diusivity (MD) increases in the right anterior thalamic radiation.

A holistic approach

Grey and white matter are thus usually studied separately and independently using either T1-w or dMRI scans respectively. However, due to the interwoven nature of the brain, the syndrome of Gilles de la Tourette is supposed to aect all the components of the circuits, namely both grey matter structures such as sub-cortical nuclei and cortical surface as well as their inter-linking white matter tracts. Associated abnormalities could alter the shape of each one of these components as well as their relative position. Following [START_REF] Dryden | Statistical shape analysis[END_REF], we dene the shape of a structure as the geometrical information invariant to the Euclidean similarity trans-1.5. A holistic approach formations. This implies that two structures have the same shape if the rst one can be transformed into the second one using only translations, scaling and rotations.

When analysing an anatomical complex, such as the CSTC circuits, abnormalities can aect also the position of a structure with respect to the other ones, namely their relative position. For instance, a dysfunction of the ux of information in the CSTC circuits can modify the location where white matter tracts integrate grey matter structures, namely the structural connectivity, as well as the spatial organisation of the sub-cortical nuclei. These anatomical alterations could be captured only analysing the CSTC circuits as a whole, that is studying both grey and white matter structures simultaneously.

As previously seen, this would imply that the information about the grey matter structures captured from T1-w scans should be combined with the information about white matter from dMRI scans. However, it is still unclear how to consistently integrate these two dierent imaging modalities into a single computational framework. Few integrative multi-modal approaches have been proposed. Two examples are [START_REF] Avants | Dementia induces correlated reductions in white matter integrity and cortical thickness: A multivariate neuroimaging study with sparse canonical correlation analysis[END_REF] and [START_REF] Savadjiev | Fusion of white and gray matter geometry: A framework for investigating brain development[END_REF] where the authors proposed to investigate the relationship between geometrical characteristics of grey and white matter using correlation analysis and mutual information respectively. These approaches can show statistically related abnormalities in grey and white matter which can help to clarify how a certain syndrome acts on the whole brain.

However, these methods independently extract only specic morphological characteristics from each scan which are then statistically related. This does not permit to take into consideration the organisation of the circuits and the conclusions that can be drawn are limited to the morphological characteristics chosen (i.e. volume, grey matter concentration, FA).

In this thesis, we follow a dierent strategy which allows us to integrate the information of dierent modalities into a single virtual representation of the brain anatomy. We propose to shift the morphological analysis from images to complexes of 3D meshes segmented from the dierent structural scans. In the case of CSTC circuits, cortical surface and sub-cortical nuclei are segmented from T1-weighted scans and they are modelled as 3D surfaces. Estimates of the trajectories of large groups of white matter neural axons are instead traced from dMRI scans using tractography algorithms and they are modelled as bundles of 3D curves called streamlines. An example of anatomical complex of 3D meshes is shown in Fig. 1.2. Tractography is the only non invasive way to model the white matter architecture and therefore it is also the only technique capable to capture the information about structural connectivity. Studying how and where neural bers integrate the grey matter structures is essential in order to understand which are the functional areas implicated in the abnormal brain development of GTS and how this syndrome modify the organisation of the CSTC circuits.

The proposed approach permits to visualise and analyse all the components of the CSTC circuits as a single complex. Our goal is to verify and, if conrmed, even localize and describe the presence of anatomical abnormalities associated to GTS within the CSTC circuits. This can be achieved by comparing two groups Chapter 1. Introduction of anatomical complexes representing the CSTC circuits of controls and patients respectively. The rst step is to dene a computational model and metric for both 3D curve and surface meshes to compare the representations of the circuits across individuals. This model should then be integrated into a statistical framework to characterise the invariants and the organisational and morphological variability of each clinical group. Eventually, we should be able to compare these statistics in order to localize and quantify the potential abnormalities aecting the GTS group.

Furthermore, the use of a statistical method also increases the robustness of the results with respect to possible imperfections due to segmentation or tractography errors. Since these errors should not be systematic, they should not inuence (or at least to a smaller extent) the estimate of the common morphological and organisational characteristics as well as their variability.

It is important to say that another possible strategy would be to use an iconicgeometric approach using both images and meshes like in [START_REF] Siless | Joint T1 and Brain Fiber Log-Demons Registration Using Currents to Model Geometry[END_REF]; Graciano Fouquier et al. (2014). However, these strategies have been used so far just to improve the accuracy of image-based registrations. Moreover, results based on images and meshes would be more dicult to visualise, and therefore also to interpret, than using only meshes as proposed in this thesis.

Mesh-based morphometry

A 3-dimensional surface mesh is composed of connected polygonal faces (usually triangles) which recreate the 3D shape of the grey matter objects. Meshes are built 1.6. Mesh-based morphometry starting from the outlines of the structures on each MRI scan, namely the segmentation results, using algorithms such as marching cubes [START_REF] Lorensen | Marching Cubes: A High Resolution 3d Surface Construction Algorithm[END_REF].

Curve meshes resulting from tractography algorithms are called streamlines, they link two areas of the dMRI scans and they are constituted of segments connecting neighbouring voxels. In this thesis we will consider only streamlines connecting two grey-matter structures of the CSTC circuits.

1.6.1 Computational models for complexes of 3D meshes Dierent computational models for both surface and curve meshes have been proposed and probably the most well known is the point distribution model [START_REF] Cootes | Active Shape Models-Their Training and Application[END_REF]. It relies on the denition of consistent point correspondences (landmarks) between meshes of dierent subjects. They can be either manually identied by an expert [START_REF] Bookstein | Landmark methods for forms without landmarks: morphometrics of group dierences in outline shape[END_REF] or automatically estimated using, for instance, dynamic particles [START_REF] Cates | Particle-Based Shape Analysis of Multi-object Complexes[END_REF] or by propagating the landmarks of an atlas through an ane [START_REF] Brett | A Method of Automated Landmark Generation for Automated 3d PDM Construction[END_REF] or non-rigid transformation [START_REF] Frangi | Automatic construction of multiple-object three-dimensional statistical shape models: application to cardiac modeling[END_REF]. Landmarks have been used mainly for 3D surface meshes since nding point-correspondences between white matter ber bundles composed of hundreds of thousands of streamlines is a dicult, if not impossible, task especially when the streamlines do not share a similar length and starting/ending anatomical locations.

Several extensions of this method have been proposed and they are based on correspondences dened using other shape parametrisations such as spherical harmonics [START_REF] Brechbühler | Parametrization of Closed Surfaces for 3-D Shape Description[END_REF], minimum-description length [START_REF] Davies | A minimum description length approach to statistical shape modeling[END_REF], skeleton-based representations [START_REF] Golland | Statistical Shape Analysis Using Fixed Topology Skeletons: Corpus Callosum Study[END_REF], m-reps [START_REF] Yushkevich | Intuitive, localized analysis of shape variability[END_REF], cubic B-splines [START_REF] Corouge | Fiber tract-oriented statistics for quantitative diusion tensor MRI analysis[END_REF] and Fourier descriptors [START_REF] Székely | Segmentation of 2-D and 3-D objects from MRI volume data using constrained elastic deformations of exible Fourier contour and surface models[END_REF][START_REF] Batchelor | Quantication of the shape of ber tracts[END_REF]. These shape descriptors have been broadly and successfully employed to analyse several brain structures in both single-and multi-object studies. However, most of them are conceived for only a particular kind of mesh, i.e. genus-zero surfaces [START_REF] Lu | Statistical Multi-Object Shape Models[END_REF][START_REF] Gorczowski | Multi-Object Analysis of Volume, Pose, and Shape Using Statistical Discrimination[END_REF] or streamlines [START_REF] O'donnell | Unbiased Groupwise Registration of White Matter Tractography[END_REF], which prevents the integration of all the components of the CSTC circuits in a unied framework.

In contrast to previous shape parametrisations, the frameworks of currents [START_REF] Glaunès | Transport par diéomorphismes de points, de mesures et de courants pour la comparaison de formes et l'anatomie numérique[END_REF][START_REF] Vaillant | Surface Matching via Currents[END_REF][START_REF] Bibliography Durrleman | Statistical models of currents for measuring the variability of anatomical curves, surfaces and their evolution[END_REF] and varifolds (Charon and Trouvé, 2013b) are not limited to a specic kind of mesh and they do not need point-correspondences or streamline-correspondences when dealing with white matter ber bundles. The framework of varifolds can be seen as the non-oriented version of the one of currents, in the sense that meshes do not need a consistent orientation among subjects. These models can thus provide a unifying framework to process all the components of the CSTC circuits. Moreover, they embed all structures in a Hilbert space which means that, for instance, the union of dierent streamlines forming a ber bundle is modelled as a sum of vectors. This simplies standard computations such as norms, distances and calculus making these frameworks two valid candidates on which ground the subsequent statistical shape analysis.

Eventually, an important characteristic of shape descriptors is their computa-Chapter 1. Introduction tional eciency, namely the computational load required to process the modelled meshes. In the case of white matter ber bundles this is fundamental since the huge number of streamlines resulting from tractography algorithms can make any analysis infeasible. To this end, dierent parsimonious representations have been proposed (see Chapter 3 for a small review) such as the one in Durrleman et al.

(2011a) based on a matching-pursuit algorithm for currents. This approximation is very concise but it has the drawback to accurately approximate only the areas of the bundle characterised by a high density of streamlines, like the central mass, almost ignoring the external smaller fascicles. Moreover, this approximation does not conserve the topology of streamlines since it can be seen as a set of disconnected points. This prevents the study of structural connectivity and therefore it can not be used to test the aforementioned hypothesis about GTS.

Statistical shape analysis of complexes of 3D meshes

Once described how to mathematically represent the structures of an anatomical complex and how to measure a distance between them, one needs to dene how to statistically analyse their morphological and organisational variability within each clinical group. Dierences between the two groups will characterise the abnormalities related to GTS.

The rst step of a statistical shape analysis is to globally align all the complexes using rigid or ane transformations in order to lter out the uninformative shape variability 1 . Transformations are applied to all the structures of the complex at the same time and the remaining dierences regard the local relative position between components and their individual shape. The alignment can be obtained using, for instance, Procrustes analysis [START_REF] Dryden | Statistical shape analysis[END_REF] which is based on the minimisation of the sum-of-squared distances between the structures with respect to the similarity transforms chosen. Distances are dened as the standard Euclidean metrics when dealing with landmarks and as geodesic distances based on appropriate Riemannian metrics when dealing with other shape descriptors such as m-reps [START_REF] Gorczowski | Multi-Object Analysis of Volume, Pose, and Shape Using Statistical Discrimination[END_REF].

Once all complexes are in the same reference frame it is possible to compute the mean conguration and the main morphological variations of the population based on a Principal Component Analysis (PCA) when working with Euclidean vector spaces (i.e. landmarks) or Principal Geodesic Analysis (PGA) when complexes lie on nonlinear Riemannian spaces [START_REF] Fletcher | Principal geodesic analysis for the study of nonlinear statistics of shape[END_REF]. Average congurations and main modes of variations of the two clinical groups can be compared in order to localize and quantify the abnormalities. Otherwise, dierences can also be visualised by deforming the average conguration of one clinical group towards the orthogonal direction to the most discriminative hyperplane in the feature space, namely the direction that best separates the two clinical groups. The feature space can be composed by the parameters describing the shape of each individual component 1 Depending on the hypothesis of the study, scaling can be an important feature.

1.6. Mesh-based morphometry and/or their relative position, depending on the hypothesis to test. This strategy is not computationally intensive and the results are easy to interpret but it has the drawback that it does not take into consideration the fact that the anatomical organisation of the circuits should not change, namely that separated structures should not intersect, when computing the average, the main morphological variations or when deforming the average conguration towards the most discriminative direction.

From a dierent perspective, rooted in the seminal book of D'Arcy Wentworth Thompson (D'Arcy [START_REF] Wentworth | On growth and form[END_REF], Grenander proposed to quantify the morphological dierences between two anatomical complexes as the amount of deformation needed to warp the rst one to the second one. Moreover, he proposed also a generative statistical model where every complex of a population is modelled as a deformation of a common anatomical complex called template complex [START_REF] Grenandner | General pattern theory -A mathematical study of regular structures[END_REF][START_REF] Miller | Mathematical textbook of deformable neuroanatomies[END_REF][START_REF] Christensen | Deformable templates using large deformation kinematics[END_REF][START_REF] Grenander | Computational Anatomy: An Emerging Discipline[END_REF][START_REF] Allassonnière | Towards a Coherent Statistical Framework for Dense Deformable Template Estimation[END_REF]. Deformations act on the entire ambient space and they modify the template complex such that it resembles to the subject complexes putting into correspondence the components of the template complex with the homologous ones of the subjects. Similarity between correspondent structures is dened by the metric of the computational models previously described. First examples used a xed template chosen a priori like the Talairach atlas [START_REF] Davatzikos | A computerized approach for morphological analysis of the corpus callosum[END_REF] or a randomly-selected subject [START_REF] Csernansky | Hippocampal morphometry in schizophrenia by high dimensional brain mapping[END_REF]. However, this kind of templates might have imperfections due to segmentation errors and they might not be centred with respect to the population under study, biasing consequently the whole analysis. Later works have thus proposed to automatically estimate the template without imposing it at the beginning [START_REF] Joshi | Landmark matching via large deformation dieomorphisms[END_REF][START_REF] Qiu | Atlas Generation for Subcortical and Ventricular Structures With Its Applications in Shape Analysis[END_REF][START_REF] Durrleman | Morphometry of anatomical shape complexes with dense deformations and sparse parameters[END_REF]. The joint estimate of the template complex and the deformations has been called: atlas construction.

In this setting, the template complex captures the common morphological characteristics of the population, and it can be thought as an average. The deformations model instead how each subject varies with respect to this average and they represent therefore the morphological variability of the population. It is important to highlight that this framework is generic and it can be applied with any kind of mesh and shape parametrisation, on condition that they are compatible with the type of deformation employed, and even with T1-w images like in [START_REF] Joshi | Landmark matching via large deformation dieomorphisms[END_REF]; Avants and Gee (2004); [START_REF] Bhatia | Consistent groupwise non-rigid registration for atlas construction[END_REF]; [START_REF] Zhang | Bayesian Estimation of Regularization and Atlas Building in Dieomorphic Image Registration[END_REF]. This avoids the use of segmentation and tractography algorithms but it prevents the analysis of structural connectivity since it is impossible to understand where white matter bers integrate grey matter structures from T1-w scans. Furthermore, in the case of 3D surfaces and 2D/3D curves (open and closed) other authors have also proposed to dene the deformations directly on the meshes rather than assuming an underlying deformation of the whole space [START_REF] Joshi | A Novel Representation for Riemannian Analysis of Elastic Curves in Rn[END_REF][START_REF] Michor | An overview of the Riemannian metrics on spaces of curves using the Hamiltonian approach[END_REF][START_REF] Younes | A metric on shape space with explicit geodesics[END_REF][START_REF] Kurtek | Parameterization-Invariant Shape Comparisons of Anatomical Surfaces[END_REF][START_REF] Bauer | A New Riemannian Setting for Surface Registration[END_REF][START_REF] Bauer | Constructing reparameterization invariant metrics on spaces of plane curves[END_REF][START_REF] Gutman | A Riemannian Framework for Intrinsic Comparison of Closed Genus-Zero Shapes[END_REF]. However, these methods have been usually designed for single-object studies and only for one specic kind of mesh. Thus, they can not be Chapter 1. Introduction employed for the analysis of the entire CSTC network.

In order to preserve the anatomical organisation of the template complex, deformations are usually dened as dieomorphisms, namely smooth bijective deformations with smooth inverse 2 . This kind of deformation prevents the structures to intersect, fold, torn or collapse to single points. This is very important in order to obtain an anatomically relevant analysis.

The group of dieomorphisms constitutes a non-linear curved manifold which is not a vector space. Statistics such as average or principal modes need to be dened in the tangent space of this manifold [START_REF] Vaillant | Statistics on diffeomorphisms via tangent space representations[END_REF][START_REF] Pennec | Intrinsic Statistics on Riemannian Manifolds: Basic Tools for Geometric Measurements[END_REF]. In the framework known as large deformation dieomorphic metric mapping (LDDMM) [START_REF] Dupuis | Variational problems on ows of dieomorphisms for image matching[END_REF][START_REF] Trouvé | Dieomorphisms Groups and Pattern Matching in Image Analysis[END_REF][START_REF] Beg | Computing Large Deformation Metric Mappings via Geodesic Flows of Dieomorphisms[END_REF], dieomorphic deformations are dened as the last ones of a ow produced by integrating a time-varying vector eld which is constrained to be smooth. Based on the metric for such a group of dieomorphisms dened in [START_REF] Trouvé | An Approach of Pattern Recognition through Innite Dimensional Group Action[END_REF], the ow that warps the template complex to every subject complex is chosen as the geodesic path between the identity transform and the last deformation of the ow (i.e. the path with the minimal cost required to make the template similar to the subject complex). In [START_REF] Miller | Geodesic Shooting for Computational Anatomy[END_REF] it is shown that the (nal) deformation of the template complex depends only on the dual of the initial time-varying vector eld, called initial momentum, which belongs to a vector space. In this space, it is therefore possible to use linear statistical analysis such as PCA (Principal Component Analysis) to study the non-linear shape space dened as the orbit of the template complex under the action of the group of dieomorphisms [START_REF] Vaillant | Statistics on diffeomorphisms via tangent space representations[END_REF]. Moreover, it oers also a computationally ecient way to compute the deformations, called geodesic shooting [START_REF] Miller | Geodesic Shooting for Computational Anatomy[END_REF]. Dierent algorithms have been proposed based on this principle such as in [START_REF] Cotter | Singular solutions, momentum maps and computational anatomy[END_REF] In this thesis, we will use the one proposed in [START_REF] Durrleman | Optimal Data-Driven Sparse Parameterization of Dieomorphisms for Population Analysis[END_REF] where the time-varying vector elds are dened with a small set of control points. The number of control points, chosen by the user, determines the number of initial momenta and it does not increase with the number of analysed structures in the complex. This makes such a model well suited for the statistical analysis of anatomical complexes like the CSTC circuits. Moreover, it has already been used with both the framework of currents (Durrleman et al., 2011a) and varifolds [START_REF] Durrleman | Morphometry of anatomical shape complexes with dense deformations and sparse parameters[END_REF]) on which we will ground our analysis in the following chapters. Nevertheless, it is important to notice that for the scope of this thesis any parametric dieomorphic deformation model could have been used (provided that it allows to compute statistics).

For instance, another powerful setting which allows to compute statistics on groups of dieomorphisms is the one proposed in [START_REF] Arsigny | A Log-Euclidean Polyane Framework for Locally Rigid or Ane Registration[END_REF] where the authors parametrised dieomorphic deformations with smooth stationary velocity elds (SVF). They presented a Log-Euclidean framework which permits to com-2 Topology preservation is assured even by a homeomorphism (that is a continuous deformation with a continuous inverse). Every dieomorphism is a homeomorphism but with an addition constraint on the regularity (smoothness) of both the transformation and its inverse.

1.6. Mesh-based morphometry pute the exponential map and the logarithm map on the manifold of the group of dieomorphisms using the scaling-and-squaring method. In this setting the statistical analysis is performed in the vector space given by the logarithms of the dieomorphisms. Several extension and applications have been proposed like in All the previous dieomorphic models provide dense deformations which vary locally across the ambient space permitting to capture the variations in the relative position between separated structures such as the sub-cortical nuclei [START_REF] Wang | Validity of large-deformation high dimensional brain mapping of the basal ganglia in adults with Tourette syndrome[END_REF]. However, using a single dieomorphism we implicitly assume that the relative position between structures in contact with each other does not change across subjects. This implies that every white matter ber tract of the CSTC circuits should link the same areas of the cortical surface and basal ganglia across the whole population. This assumption is in contradiction with the aforementioned hypothesis about GTS and it precludes the study of abnormalities in structural connectivity which could be caused by an abnormal brain development.

GTS should be studied with an holistic approach where morphometry and structural connectivity analysis are unied into a single computational framework. However, structural connectivity is usually studied independently of the integrated grey matter structures without considering neither their morphology nor their organisation.

Structural connectivity analysis

Structural connectivity refers to the density of connections resulting from a tractography method and integrating grey matter structures.

The most popular strategy to study dierences in structural connectivity between two clinical groups is based on the subdivision of the cortical surface and sometimes even sub-cortical nuclei in parcels which are reproducible across the subjects of both groups. This step is usually based on the registration of all subject T1-w scans to a common atlas (i.e. the MNI-space, Montreal Neurological Institute) which is already parcelled. After that, the spatial location of the parcels in the T1-w scan is transferred to the DWI scan of the same subject through registration. Then, it is used a tractography algorithm on the parcelled DWI scans to trace streamlines from specic ROI (Region of Interest) or throughout the whole brain. Eventually, one counts the number of streamlines connecting two ROIs, two parcels or a determined ensemble of parcels of the cortical surface (i.e. motor area) to a sub-cortical nucleus. The density of streamlines quanties the structural connectivity. Multi-Chapter 1. Introduction variable or uni-variate statistics can then be used to nd the most discriminating parcels between the two groups. Discordant results exist about GTS. For instance in [START_REF] Makki | Altered Fronto-Striato-Thalamic Connectivity in Children with Tourette Syndrome Assessed with Diusion Tensor MRI and Probabilistic Fiber Tracking[END_REF]; [START_REF] Cheng | Altered intrahemispheric structural connectivity in Gilles de la Tourette syndrome[END_REF] the authors have found that in patients both caudate and thalamus have a reduced connectivity with the anterior dorsolateral frontal cortex and that the frontal cortical areas have a lower connectivity with the basal ganglia. On the other hand, in [START_REF] Worbe | Altered structural connectivity of cortico-striatopallido-thalamic networks in Gilles de la Tourette syndrome[END_REF] the authors pointed out that both the motor circuit (primary motor and sensory cortical areas) and the associative circuit (lateral orbitofrontal cortex) presented an abnormally enhanced structural connectivity with the striatum and thalamus in patients subject to GTS.

The discrepancy between these studies could be due to methodological dierences or to the small size of the data samples in the rst two studies (15 and 18 subjects respectively with respect to 77 in [START_REF] Worbe | Altered structural connectivity of cortico-striatopallido-thalamic networks in Gilles de la Tourette syndrome[END_REF]). Nevertheless, this also indicates that there is a need of new, rigorous and robust statistical methods to brighten the GTS mechanism aecting structural connectivity.

The eld of complex networks [START_REF] Bullmore | Complex brain networks: graph theoretical analysis of structural and functional systems[END_REF][START_REF] Rubinov | Complex network measures of brain connectivity: Uses and interpretations[END_REF], based on graph theory, proposes innovative solutions to study structural connectivity as a network where the nodes are the parcels and the edges are constituted by the number of streamlines (or other indices) linking two nodes. Dierent measures characterising a network exist such as the degree distribution, small-world topology, highly connected hubs and modularity. They usually have a physical meaning which makes them easily interpretable. There are already some works about GTS which make use of these graph measures [START_REF] Worbe | Functional immaturity of cortico-basal ganglia networks in Gilles de la Tourette syndrome[END_REF] but most of them are about functional networks and not about structural connectivity networks. Moreover, to the best of our knowledge, there are no methods based on the complex network theory which integrate morphology and structural connectivity analysis.

Contributions and Manuscript Overview

The purpose of this thesis is to provide a methodology which is able to discover and characterise morphological and organisational abnormalities in anatomical complexes composed of both grey and white matter structures. To this end, the complexes must be studied as a whole and not as a set of independent components.

The proposed method should be able to discover both local and distributed shape alterations across the whole anatomical complex. Moreover, it should unveil also pathological changes altering the relative position of both separated components (i.e. sub-cortical nuclei) as well as connected structures (i.e. structural connectivity). All these abnormalities should be quantied, accurately localised and easily visualised in order to be useful in the clinical domain. The proposed method has been driven by the will to test the aforementioned hypothesis about GTS. Nevertheless, it could also be applied to the study of other neurodevelopmental syndromes which are thought to be associated with alterations of the CSTC circuits or other neural networks.

Contributions and Manuscript Overview

We propose to ground our method on the generative model of the previously described Grenander's atlas construction procedure. This allows us to estimate the common features of a population of anatomical complexes, represented by the template complex, and at the same time the morphological variability described by the set of dieomorphisms, one for each subject. These characteristics can then be compared between a group of healthy subjects and one of GTS patients in order to localize the possible abnormalities due to the syndrome. This strategy has been mainly employed for studying single structures or complexes composed of separated structures. Here we propose to extend this method to anatomical complexes composed of grey matter structures interconnected by white matter tracts. This kind of anatomical complexes can not be studied with previous deterministic, single-dieomorphic atlas procedures based on computationally expensive or topology-neglecting white matter shape models.

To this end, we propose to embed the atlas construction in a Bayesian framework which permits to automatically estimate important balancing parameters, one for each structure of the complex, which were xed by the user in previous methods and which highly inuenced the results. Moreover, this statistical setting also permits to directly estimate a consistent and well-conditioned covariance matrix, which is fundamental for the statistical analysis of the morphological and organisational variability. In order to capture the variations in the relative position of both separated and attached structures, making thus possible the analysis of structural connectivity, we propose also an innovative deformation setting called double diffeomorphism. It is a cascade of two dieomorphisms which can model the variation in relative position between white and grey matter structures and at the same time correctly warp the template complex to every subject complex. Eventually, to reduce the computational resources required to work with white matter ber bundles, we also propose to approximate them with a parsimonius representation based on weighted prototypes. This makes feasible the computation of a Bayesian double dieomorphic atlas. These contributions are thoroughly explained in the following chapters. It follows a brief description of each one of them.

1. Chapter 2: The estimate of the atlas is usually based on a deterministic approach where a gradient descent scheme is used to minimize a cost function constituted of M data-terms, one for each structure of the complex, and a regularity term on the dieomorphisms. The inuence of the data-terms is balanced by M trade-o parameters which are xed at the beginning. The choice of these values is crucial since they inuence the results. A value too small might weight too much the data term leading to a situation of overtting. On the contrary, a value too large might penalize the deformations reducing the accuracy of the matching. Moreover, the number of deformation parameters (initial momenta) is usually much higher than the number of subjects which makes the atlas construction a high-dimensional low samples size (HDLSS) problem. The study of the morphological variability of a population Chapter 1. Introduction is often based on a Principal Component Analysis (PCA) which needs the computation of the covariance matrix of the deformation parameters. In a HDLSS setting one needs to regularize this estimate at the end of the atlas construction making therefore the computed covariance matrix inconsistent with the optimal deformations resulting from the atlas procedure. To this end, we propose to embed the atlas construction in a Bayesian framework taking advantage of the prior distributions introduced in [START_REF] Allassonnière | Towards a Coherent Statistical Framework for Dense Deformable Template Estimation[END_REF]. This allows us to estimate both the trade-o parameters and a wellconditioned covariance matrix of the deformation parameters which can be directly employed in statistical analysis such as PCA. Furthermore, we propose to model all meshes as Gaussian random varifolds. Since the space of varifolds is of innite dimension, we dene a grid for every structure, similarly to [START_REF] Bibliography Durrleman | Statistical models of currents for measuring the variability of anatomical curves, surfaces and their evolution[END_REF], on which both shapes and template are projected to. This determines a nite-dimensional approximation space on which it is possible to dene probability density functions. Eventually, we extend also the proposed Bayesian framework to a multi-population atlas construction.

We take advantage of the clinical diagnosis of the subjects to divide them into dierent groups (i.e. controls and patients). The groups share the same template and their deformation parameters are modelled with dierent distributions which are dened on the same space. This allows us to compare them and to test hypothesis of equality between the moments of their distributions such as their covariance matrices.

2. Chapter 3: Fiber bundles stemming from tractography algorithms contain many streamlines, from few thousands up to millions. They require therefore a great amount of computer memory and computational resources to be stored, visualised and processed. This can limit their use in computationallyintensive processes such as the atlas construction. To this end, we propose to approximate a ber bundle with a parsimonious representation of weighted prototypes. Prototypes are chosen among the streamlines and they represent groups of similar streamlines. Their weight is related to the number of streamlines approximated. Both streamlines and prototypes are modelled with a new extension of the framework of currents, which we have called weighted currents. This computational model inherits all the qualities of usual currents, such that it does not need point-to-point correspondences, and it augments its geometry-based metric with information about structural connectivity in the spirit of functional currents (Charon and Trouvé, 2013a). Two streamlines are considered similar if their endpoints are close to each other and if their pathways are similar. Moreover, the space of weighted currents is also a vector space with a closed-form metric. This permits to easily compute the approximation error and to select the prototypes based on the minimisation of this error. We propose an iterative algorithm which approximates independently and simultaneously all the fascicles of the bundle. This allows us to approximate the whole bundle in a fast and accurate way. We evaluate our algorithm 1.7. Contributions and Manuscript Overview on bundles resulting from both deterministic and probabilistic tractography algorithms. We show that this approximations use on average only 2% of the original streamlines as prototypes. This drastically reduces the computational burden required for any processing such as atlas construction while preserving both the geometry and the structural connectivity of the original bundles. The latter quality is fundamental for the scope of this thesis.

3. Chapter 4: Single-dieomorphic deformations can correctly put into correspondence anatomical complexes composed of separated structures. Since they preserve topology, they do not alter the organisation of the anatomical complexes preventing structures to fold, intersect or shear. However, this also implies that the relative position between structures in contact with each other should not change across subjects. This means, for instance, that white matter ber bundles should link the same functional areas of the cortical surface and basal ganglia across the whole population. This assumption precludes the study of changes in structural connectivity and therefore it also impedes to test the aforementioned hypothesis about GTS. We propose to tackle this problem with a new deformation setting called double dieomorphism. The template complex is warped towards every shape complex of the population using a cascade of two dieomorphisms. The rst one acts only on the white matter of the template complex, keeping xed the grey matter. It can be seen as a relative change of coordinates within the gray matter considered as a xed reference frame. During this transformation the ber bundles slide on the xed gray matter structures, capturing the variation in structural connectivity. These variations can be compared among subjects since they are all computed with respect to the same reference frame, namely the xed gray matter of the template complex. The second dieomorphism acts on the whole template complex, namely on both the resulting white matter and on the gray matter. This is a global change of coordinates which has the role to put into correspondence both gray and white matter structures of the template complex with the homologous ones of the subject's shape complex. The rst dieomorphism re-positions the white matter tracts within the gray matter such that all the components of the template complex can be correctly put into correspondence by the second dieomorphism. The two dieomorphisms are optimised together minimising a single cost function and using a gradient descent scheme. This new deformation setting is combined with the Bayesian framework presented in Chapter 2 into a Bayesian double dieomorphic atlas construction. We use this framework to study the morphological and organisational variations of the cortico-putamen circuit in both controls and GTS patients. By comparing these results we are able to visualise and quantify the dierences between the two clinical groups both in terms of shape and structural connectivity. The feasibility of this study is due to the parsimonious representation for white matter ber bundles proposed in Chapter 3.

Eventually, we also show that a double dieomorphic approach leads to better Chapter 1. Introduction classication scores than using the single dieomorphic deformation presented in Chapter 2. We use 49 GTS patients and 27 healthy controls with a Linear Discriminant Analysis (LDA) and a leave-one-out cross-validation strategy. In this chapter, we focus on a general morphometric approach called atlas construction. It consists of estimating an average shape complex of a population called template complex and the deformations of the embedding space which warp the template complex to the shape complexes of every subject. Deformations put into correspondence the template with every shape complex. The template captures the common morphological characteristics of the population and the deformations capture the variations in shape across the population. Every deformation is a single 3D dieomorphism (a smooth deformation with smooth inverse) which transforms the entire 3D space. All meshes are therefore deformed together whatever their number or kind. It permits to preserve the spatial organization of the structures since neither intersection, folding or shearing may occur. This allows us to study a human organ as a whole and not as a set of independent components, making thus possible to take into account also the information about the relative position between separated structures of the template complex.

The deformations transform the template complex so that each one of its structures superimposes with the homologous one of the subject shape complex. To this end, one needs to dene a similarity measure for every structure. Dierent metrics have been proposed for surfaces and curves. Some of them assume that it is possible to nd correspondences between points of homologous structures, called landmarks. This is usually hard (or even impossible) especially when working with groups of curves. This is why we have opted for a metric based on the framework of varifold Here we propose to embed a dieomorphic multi-object atlas construction into a Bayesian framework. Shape complexes can be composed of curves and surfaces and they are both modelled as varifolds. The noise of each structure is assumed to be a random Gaussian varifold which belongs to an innite dimensional space. We dene a nite-dimensional approximation space where the noise is projected to and where it is possible to dene probability distributions. Furthermore, we use a prior distribution on the variance of the noise as in Allassonnière et al. (2007) adapted for shape complexes modelled as varifolds. This permits to automatically estimate the balancing weights in one single atlas construction. Moreover, we assume that also the deformation parameters follow a Gaussian distribution. We use a prior distribution on their covariance matrix similarly to [START_REF] Allassonnière | Towards a Coherent Statistical Framework for Dense Deformable Template Estimation[END_REF] which allows to estimate it automatically even with a small sample size. Furthermore, the estimated covariance matrix is always well-conditioned whatever the number of subjects and it can be directly employed to compute statistics about the population under study.

In addition to that, we also take advantage of the proposed Bayesian framework to dene a multi-population atlas construction. We use the information given by the clinical diagnosis of the subjects to subdivide the data-set in dierent groups.

Their deformation parameters are then modelled as dierent Gaussian distributions.

Since all sets of deformation parameters are dened on the same space, we can quantitatively compare them and the moments of their distributions such as the covariance matrices. The proposed Bayesian framework permits to perform this multi-population analysis at no additional cost, namely without tuning any new balancing weight.

The chapter is organised as follows. In Section 2.2, we rst formulate the proposed Bayesian framework in a rather general way, showing that it could be employed with dierent shape models and deformation settings. We then present how to model multi-object shape complexes with varifolds and how to dene a varifold random variable. Afterwards, we show the dieomorphic framework used in the 2.2. Methods experiments and how to integrate it in our statistical setting. Then, we present the gradient-descent scheme used to optimise the cost function and how to initialise the template of both surfaces and curves. In the last part of Section 2.2, we describe how to compute a multi-population atlas construction. Eventually, we test the robustness of our algorithm with respect to the hyperparameters of the prior distributions and we show also how the proposed method can be employed to assess morphological dierences between two groups of subjects.

Methods

Bayesian atlas construction

We dene the mesh of structure j belonging to subject i as S ij . The total number of subjects is N and the number of analysed structures is M . Every subject shape complex S i , dened for the moment in a generic way as the ensemble of all the meshes S ij of subject i, is modelled as a deformed template complex φ i (T ) plus a residual ε i . Both T and ε i are also dened as the ensembles of the templates T j and residuals ε ij . This formulation is also known as forward model where we assume that all elements belong to an algebraic structure where addition is dened

S i = φ i (T ) + ε i (2.1)
The deformation φ i , proper to subject i, can belong to any dieomorphic framework present in the literature. The only requirement, for the scope of this chapter, is that it deforms the entire ambient space, namely all the structures of the template simultaneously. The goal of the atlas construction is to estimate the template and its morphological variations within the population of shape complexes. The variations are described by the ensemble of deformations {φ i } and each one of them is parametrised by a set of parameters α i ∈ R q . We assume that these parameters follow a Gaussian distribution with zero mean and covariance matrix Γ α :

α i ∼ N (0, Γ α ) p(α i |Γ α ) ∝ 1 |Γ α | 1/2 exp - 1 2 α T i Γ -1 α α i (2.2)
The mean is set to zero because the template is supposed to represent the average of the shape complexes. In other words, using the terminology of mixed models, the deformation parameters α i are the random eects associated to the xed eect T which is therefore the expected value of the shape complexes.

From Eq.2.2 we can notice that the distribution of each set of deformation parameters α i is completely described by Γ α and therefore we can rephrase our goal as estimating the template T and the covariance matrix Γ α knowing the shape complexes of the population {S i }, or more formally maximizing their joint posterior Chapter 2. A Bayesian framework for multi-object atlas construction distribution:

{T * , Γ * α } = arg max T ,Γα p(T , Γ α |{S i }) (2.3)
This maximization is constrained by the fact that the deformed template φ i (T ) should resemble to the shape complex S i or, in other words, that the residual ε i should be small. This is a common problem in statistical learning and it is usually tackled by assuming that the residual follows a Gaussian distribution centred at 0. The maximization of its likelihood is equivalent to minimize its squared norm. We will start by modelling the structures of the shape complexes S i and of the template complex T with landmarks. Every structure j has a number of landmarks equal to Λ j and the norm of the dierence between two meshes of the same structure is dened as the square root of the sum of squared dierences between pair of landmarks (L 2 -norm, || • || 2 ). The likelihood of the residuals modelled as landmarks is thus dened as:

ε ij ∼ N (0, σ 2 j Id Λ j ) p(ε ij |σ 2 j ) ∝ 1 |σ 2 j | Λ j /2 exp - 1 2σ 2 j ||S ij -φ i (T j )|| 2 2 (2.4)
In section 2.2.2 we will make clear how to adapt this equation to varifolds by changing only the norm || • || and the denition of the parameter Λ j . The variance of the noise of structure j depends only on the parameter σ 2 j which is structure-dependent and it can also be modelled as a random variable. Moreover, it is possible to notice that also the shape S ij follows a Gaussian distribution: S ij ∼ N (φ i (T j ), σ 2 j Id Λ j ). Assuming independence between all random variables and considering α i as an unobserved nuisance variable we can rewrite Eq.2.3 as:

{T * , Γ * α , σ 2 * j } = arg max T ,Γα,σ 2 j N i M j p(T j , Γ α , σ 2 j , α i , S ij )dα i (2.5)
The maximization of Eq.2.5 is not tractable analytically. A possible solution is to employ an iterative algorithm like the EM (Expectation Maximization) and to approximate the conditional distribution of the E step with a Dirac distribution at its mode. Let Θ = {Γ α , T , {σ 2 j }} be the parameters of interest and Z = {{α i }} the unobserved nuisance variables, it results:

{Θ * , Z * } = arg min Θ,Z -log [p({S i }/Z, Θ)] -log [p(Z/Θ)] -log [p(Θ)] (2.6)
This algorithm is an approximation of the MAP estimator which does not converge when applied to a dataset with low signal-to-noise ratio (SNR), as demonstrated in [START_REF] Allassonnière | Towards a Coherent Statistical Framework for Dense Deformable Template Estimation[END_REF]. In order to ensure the convergence of the iterative estimation scheme one may use sampling algorithms like MCMC (Markov

Methods

Chain Monte Carlo) [START_REF] Allassonnière | Construction of Bayesian deformable models via stochastic approximation algorithm: a convergence study[END_REF]. Unfortunately, these methods require a great computational load and execution time.

Not using priors for σ 2 j and Γ α is equivalent to a ML estimation and it can produce degenerate estimates. Instead, as demonstrated in Allassonnière et al.

(2007), the introduction of inverse Wishart distributions makes possible to obtain good estimates even when the number of subjects is small. They regularize the estimates of σ 2 j and Γ α . Their probability density functions are:

σ 2 j ∼ W -1 (P j , w j ) p(σ 2 j ; P j , w j ) ∝ (σ 2 j ) -w j 2 exp - 1 2
w j P j σ 2 j (2.7) Γ α ∼ W -1 (P α , w α ) p(Γ α ; P α , w α ) ∝ |Γ α | -wα 2 exp - 1 2 w α T r(P T α Γ -1 α ) (2.8)
The scalars w j , P j and w α are strictly positive and P α is a positive symmetric matrix. Using these priors for σ 2 j and Γ α and assuming that the template T has a non-informative prior distribution, Eq.2.6 becomes:

M j=1 N i=1 1 2σ 2 j ||S ij -φ i (T j )|| 2 2 + P j w j N + M j=1 1 2 (w j + Λ j N ) log(σ 2 j )+ 1 2 N i=1 (α i ) T (Γ α ) -1 α i + 1 2 (w α + N ) log(|Γ α |) + w α 2 tr((Γ α ) -1 P α ) (2.9)
The framed terms refer respectively to the data-terms and to the regularity of the deformations. If we set w j = 0, their sum is the cost function minimised in standard atlas construction procedures Avants and Gee (2004); Durrleman et al.

(2014) where it is not employed a statistical setting. This minimisation is often referred to as the computation of the Fréchet mean. Both σ 2 j and Γ α are xed by the user and it can be noticed that σ 2 j weights the contribution of structure j as mentioned in the Introduction.

The proposed statistical framework is general and it could be applied to any shape model and dieomorphic setting provided that one denes a probability density function for the noise and a parametric deformation model. Modelling shapes with landmarks simplies the computation of the metric and of the gradient and it eases also the denition of random variables. Unfortunately, the choice of corresponding points between homologous structures of dierent subjects is not always an easy task and in some cases, like for ber bundles in the brain, it is almost impossible. It is better to opt for a correspondence-free metric like the one based on the framework of varifolds Charon and Trouvé (2013b).

Chapter 2. A Bayesian framework for multi-object atlas construction

Varifolds

The framework of varifolds is an extension of the one of currents Vaillant and Let X and Y be two curves (oriented or unoriented rectiable curves), they can be modelled as varifolds based on:

V X (ω) = X ω(x, ← -→ β(x))|β(x)| 2 dx V Y (ω) = Y ω(y, ← -→ γ(y))|γ(y)| 2 dy (2.10)
where V X (ω) and V Y (ω) are scalars, β(p) and γ(y) are the tangents of X and Y at the point x and y respectively and ← → β p is the unoriented unit vector with the same direction of β p (respectively ← → γ q for γ q ). The set of unoriented unit vectors ← → β ∈ ← → S can be formally dened as the quotient of the unit sphere in R 3 by the two elements group {±Id 3 }. Instead, the test eld ω ∈ W is a function of both the points in the ambient space R 3 and the unoriented unit tangent vectors in ← → S . The curves X and Y can be seen as polygonal lines of P and Q segments respectively. Every segment of X is completely described by its center point x p ∈ R 3 and tangent vector β p ∈ R 3 centred at x p (respectively y q ∈ R 3 and γ q ∈ R 3 for Y ). Let c p be the length of the tangent vector β p (respectively d q for γ q ), one can approximate Eq.2.10 with: V X (ω)

≈ P p=1 ω(x p , ← → β p )c p and V Y (ω) ≈ Q q=1 ω(y q , ← → γ q )d q .
It is important to notice that, even if this approximation seems reasonable, it is still an open question how to ensure its convergence when the number of segments tends to innite, i.e. with a more accurate sampling of the curve [START_REF] Charlier | The fshape framework for the variability analysis of functional shapes[END_REF].

Furthermore, Eq.2.10 shows that the space of varifolds W * is a linear functional from W to R. This makes W * a vector space and therefore the union of meshes is equal to a sum in the space of varifolds. But, contrary to currents, inverting the orientation of a mesh does not change the varifold representation. One chooses W to be a separable RKHS on R 3 x ← → S in order to have an explicit denition of distance between two meshes modelled as varifolds. More precisely, since W is a product space, one denes its kernel K W as the tensor product between two kernels k x , Gaussian, and k β , Cauchy-Binet, dened on R 3 and on ← → S respectively. Thus,

Methods

the inner product between X and Y modelled as varidolds is dened as:

< V X , V Y > W * = P p=1 Q q=1 k x (x p , y q )k β ( ← → β p , ← → γ q )c p d q = P p=1 Q q=1 exp -||x p -y q || 2 λ 2 W β T p γ q c p d q 2 c p d q (2.11)
The distance between V X and V Y is therefore:

||V X -V Y || 2 W * = < V X , V X > W * + < V Y , V Y > W * -2 < V X , V Y > W * .
This metric is completely parametrized by the standard deviation λ W of the Gaussian kernel k X which is xed by the user. This framework can be easily extended also to surfaces by replacing x p (y q ) and β p (γ q ) with the centres and normals of the faces of the surface, as shown in Chapter 4 and in Charon and Trouvé (2013b). The computation of the gradient with respect to the points of the mesh can be found in Charon and Trouvé (2013b) Varifolds random variable In Eq.2.4 we dened the residual ε as a multivariate Gaussian random variable on the nite-dimensional Euclidean space of landmarks.

Here we extend this denition to the framework of varifolds following the same line of reasoning presented in [START_REF] Bibliography Durrleman | Statistical models of currents for measuring the variability of anatomical curves, surfaces and their evolution[END_REF] for the framework of currents. We dene a random Gaussian varifold as a linear map between every vector eld ω ∈ W to a real random Gaussian variable G(ω) such that, given two vector elds ω 1 and ω 2 , E[G(ω 1 )] = 0 and E[G(ω 1 ), G(ω 2 )] =< ω 1 , ω 2 > W . This shows that the kernel K W of the space W completely denes the covariance matrix of the Gaussian varifold.

However, since it is innite-dimensional, it has no probability density function and therefore we can not simulate instances of it. To tackle this problem, we dene for each structure j a nite-dimensional space W * Λ j on which we project the template T j and all the shapes {S ij } i=1...N modelled as varifolds. This nite dimensional space W * Λ j is dened as the span of a set of delta Dirac varifolds: Span{δ (xu, ← → β k ) } where both the points {x u } and the unoriented unit vectors { ← → β k } are constrained to belong to two predened grids, respectively Υ x and Υ β . The rst one is a linearly spaced grid in the ambient space and Υ β is a regular sampling of the half unit sphere in R 3 . The number of points of the two grids is respectively Λ j

x and Λ β and Λ j is their product where the index j refers to the j-th structure.

For every nite-dimensional space W * Λ j , it is possible to dene a block matrix

K W,Λ j whose blocks are the RKHS kernel K W,Λ j (x u , ← → β k ), (•, •) u=1,...,Λ j x k=1,...,Λ β between every possible combination of the couples {x u , ← → β k }. The matrix K W,Λ j has a di- mension of [Λ j xΛ j ]. To project a delta Dirac varifold δ (y, ← → α ) |α| onto W * Λ j , using
for instance the closest neighbour projection, we look for the closest point to y of the grid Υ x and for the closest direction to ← → α among the ones given by Υ β and Chapter 2. A Bayesian framework for multi-object atlas construction we assign the scalar |α| to that particular couple of grid points. The projection is therefore completely parametrised by the vector {|α|} of size Λ j . Using this scheme, the squared norm of a projected varifold ← → β q } respectively. Scaling the kernel K W,Λ j with a scalar σ 2 j , we can dene the likelihood of the residuals modelled as varifolds in the following way:

|| L l=1 δ (y l , ← → α l |α l |)|| 2 W * Λj is equal to u∈Υx p∈Υx k∈Υ β q∈Υ β c uk K W,Λ j (x u , ← → β k ), (x p , ← → β q ) c
p(ε ij |σ 2 j ) ∝ 1 |σ 2 j | Λ j /2 exp - 1 2σ 2 j ||S ij -φ i (T j )|| 2 W * Λj (2.12)
where it is important to notice that K W,Λ j is a constant matrix that is computed only once at the beginning of the atlas construction. This equation has the same formulation of the likelihood of the residuals modelled as landmarks in Eq.2.4. The only two dierences are the denitions of the norm and of the parameter Λ j . These are also the only changes to apply to the cost function in Eq.2.9 when using the framework of varifolds instead than landmarks.

In practice, due to the important computational load and time required to cal-

culate || • || 2 W * Λj , we choose to compute directly ||(•)|| 2 W * using a fast GPU (CUDA)
implementation. Moreover, we dene Υ j

x as a regular linearly-spaced grid containing both shapes and template of structure j and where the distance between every couple of points is given by λ W . The resulting number of points Λ j

x is used to compute Λ j =Λ j x . In this way Λ j depends only on λ W which is xed by the user and it is related to the distance at which two shapes can be compared and at the same time to the degree of detail of the anatomical features one wants to detect. In Sec.2.3.3

we show that the number of grid points does not aect the results above a certain threshold and that the proposed choice seems a good one.

Dieomorphic Transformations

We dene here how to compute the dieomorphic deformations of the template complex. Our approach relies on the Large Deformation Dieomorphic Metric Mapping (LDDMM) framework based on the control point formulation presented in Durrleman et al. (2011b). Template transformations are built by integrating a time-varying vector eld v t (x) over t ∈ [0, 1] where v t (x) represents the instantaneous velocity of every point x belonging to the ambient space at time t. Calling φ t (x) the position of a point at time t which was located in x at time t = 0, its evolution is given by: ∂φt(x) ∂t = v t (φ t (x)) with φ 0 (x) = x, which produces a ow of deformations {φ t } t∈[0,1] . Furthermore, v t belongs to a Reproducible Kernel Hilbert Space (RKHS) with Gaussian kernel K V . In [START_REF] Glaunès | Transport par diéomorphismes de points, de mesures et de courants pour la comparaison de formes et l'anatomie numérique[END_REF] it is shown that if such a vector eld is square integrable in the interval [0,1] then every deformation of the ow {φ t } t∈[0,1] is a dieomorphism. Here we dene v t using a dynamical system of C p control points 2.2. Methods c={c k ∈ R 3 } scattered in the ambient space and a set of time-varying vectors called momenta α={α k ∈ R 3 } linked to each control point. This means that the velocity of every point of the space is given by:

ẋ(t) = v t (φ t (x)) = Cp p=1 K V (x(t), c p (t))α p (t) φ 0 (x) = x(0) = x (2.13)
where φ t (x)=x(t). With this denition every φ t is a dieomorphism provided that all α(t) are square integrable for every t. The evolution of a point x depends only on its initial position x(0) = x 0 and on the evolution of the system L(t) = {c(t), α(t)}. In an atlas construction, the deformations of the template complex are the last ones (φ 1 ) of a ow of dieomorphisms. Among all the possible paths connecting φ 0 to φ 1 we use the geodesic one, which means the one that minimizes the total kinetic energy along the path:

1 0 ||v t || 2 V dt = 1 0 Cp k=1 Cp p=1 α k (t) T K V (c k (t), c p (t))α p (t)dt (2.14)
It has been demonstrated in [START_REF] Durrleman | Morphometry of anatomical shape complexes with dense deformations and sparse parameters[END_REF] that the extremal paths are such that the system L(t) satises:

ċk (t) = Cp p=1 K(c k (t), c p (t))α p (t) αk (t) = -Cp p=1 α k (t) T α p (t)∇ 1 K(c k (t), c p (t)) (2.15)
Denoting L(0)=L 0 ={c(0), α(0)} the initial condition of the system, the previous set of ODEs can be rewritten as: L(t) = F [L(t)] with L(0) = L 0 . This equation shows that the entire ow of dieomorphisms, and more precisely the last one, is completely determined by the initial state of the system L 0 . Thus, given c(0) and α(0), one rst integrates Eq.2.15 obtaining the values of c(t) and α(t) at any time t. Then, every point x in the ambient space is deformed by integrating Eq. 2.13.

In geometrical term, the value L 0 denes the initial velocity of the geodesic path in the tangent space of our group of dieomorphisms at identity. The parameters {c(0), α(0)} can be used to perform the so-called tangent-space statistics Vaillant et al. (2004).

Optimization procedure

This dieomorphic setting can be easily included into the proposed Bayesian framework. We choose to employ a unique set of control points c for the whole population and N subject-specic sets of momenta {α i } i=1,...,N . The set of initial control points and momenta {c(0), {α i (0)} i=1,...,N } = {c 0 , {α i0 } i=1,...,N } represents the deformation parameters to warp the template complex towards the shape complex of subject Chapter 2. A Bayesian framework for multi-object atlas construction i. Thus, the hidden variables in Eq.2.6 become: Z = {{α i0 } i=1,...,N , c 0 }. Assuming that all random variables are independent and that the control points c 0 have a noninformative prior distribution, the expression of the cost function in Eq.2.9 does not change except for the second framed term which becomes

1 2 N i=1 (α i0 ) T (Γ α ) -1 α i0 .
The parameters T , {α i 0 } and c 0 are minimised using an accelerated version of the line search gradient descent method based on the Nesterov's scheme Nesterov (1983). Instead, the use of conjugate priors makes possible to compute the optimal values for {σ 2 j } and Γ α in a closed form:

Γα = N i=1 (α i0 )(α i0 ) T + w α P T α (w α + N ) σ2 j = N i=1 ||S ij -φ i (T j )|| 2 W * Λj + w j P j (w j + N Λ j ) (2.16)
The rst parameter Γα is equal to a weighted sum between the sample covariance matrix of the initial momenta α i0 and the prior P α . A good choice for the prior seems to be:

P α = K -1
V , where K V is a block matrix whose blocks are 3D Gaussian kernels between two dierent control points. It is exactly the kernel of the RKHS to which belongs the velocity eld v t building the dieomorphism φ 1 . If the number of subjects N is small with respect to w α then Γα ∝ K -1

V and this means that the deformation regularity part in Eq.2.9 becomes N i=1 (α i0 ) T K V α i0 , which is the sum of the geodesic distances between the template complex and all its transformations. This is exactly the deformation regularity term used in previous atlas constructions not based on a statistical setting [START_REF] Durrleman | Morphometry of anatomical shape complexes with dense deformations and sparse parameters[END_REF]. Moreover, the use of this prior makes possible the inversion of Γ α even when the number of subjects N is smaller than the number of deformation parameters.

The second parameter σ2 j is equal to a weighted sum between the data-term of the j-th structure and the prior. This term balances the contribution of the dierent structures to the cost function (Eq.2.9) and at the same time it is a trade-o between the data-terms and the the other parts of the cost function. The automatic estimate of this parameter is rather important as shown in Sec.2.

The use of an inverse

Wishart prior helps to avoid over-tting. In fact, using a maximum likelihood estimator, the minimization process might focus only on a structure k reducing its residuals almost to zero and consequently log(σ 2 k ) →-∞. This could also result in almost ignoring all the other structures. The role played by the prior is to impose a minimum value to σ 2 k in order to avoid such a situation.

The gradients of the cost function E (Eq.2.9) with respect to T ,{α i 0 } and c 0 are respectively equal to:

Methods

∇ T k E = N i=1 1 2σ 2 k ∇ T k [D ik ] ∇ α s0 E = M j=1 1 2σ 2 j ∇ α s0 [D sj ] + Γ -1 α α s0
(2.17)

∇ c 0 E = N i=1 M j=1 1 2σ 2 j ∇ c 0 [D ij ] + w α 2 ∇ c 0 tr((Γ α ) -1 P α )
where

D ij =||S ij -φ i (T j ))|| 2 W *

Λj

. The dierentiation of the prior P α with respect to c 0 is not taken into account since its norm is negligible with respect to the one of the data term.

In order to compute the three derivatives of the data-terms D ij , we need rst to compute the deformed template complex φ i (T ) where the deformation φ i is the nal dieomorphism of the ow generated by integrating forward in time Eq.2.15 and then Eq.2.13 which can be rewritten in matrix notation as:

Ṫ i (t)=G[T i (t), L i (t)]
with T i (0)=T (0)=T . After that, we can compute the gradient of D ij with respect to the deformed template complex T i (1) = φ i (T ) and then bring back this information to t = 0 in order to update the initial template points T , control points c 0 and momenta α i0 . This is done by integrating backward in time from t=1 to t=0 a set of linearised ODEs called adjoint equations like in [START_REF] Durrleman | Morphometry of anatomical shape complexes with dense deformations and sparse parameters[END_REF]. It results:

∇ T E = N i=1 θ i (0) ∇ α i0 E = ξ α i (0) + Γ -1 α α i0 (2.18) ∇ c 0 E = N i=1 ξ c i (0)
where the auxiliary variables ξ i (t) = {ξ α i (t), ξ c i (t)} and θ i (t) satisfy the linearised ODEs: Chapter 2. A Bayesian framework for multi-object atlas construction

θi (t) = -(∂ T i G[T i (t), L i (t)]) T θ i (t) θ i (1) = 1 2σ 2 ∇ T i (1) [D i ] (2.19) ξi (t) = -(∂ L i G[T i (t), L i (t)]) T θ i (t) + d L i F [L i (t)] T ξ i (t) ξ i (1) = 0

Template initialisation

Since we use a gradient descent scheme to minimize the cost function, we need an initial template. This should not be the mesh of a subject otherwise it would inherit all the imperfections due to segmentation errors and data noise, biasing consequently the atlas construction.

For surfaces, we use a regular sphere which is rst centred in the middle of the population and then it is scaled to an ellipsoid by using the three main modes of variability of the surfaces of the population considered as a points cloud.

For bundles of curves, we propose instead to initialise the template as the most relevant curves in the population.

We rst gather the curves F of every subject bundle into a raw initial template T = ∪ P h=1 F h , where P is the total number of curves. Then we select the set of H curves which resembles the most to the raw template T . Calling {F } k one of the possible P H

combinations of H curves out of P , it results:

arg min k ||T -{F } k || 2 W * k = 1, ..., P H (2.21)
which means that we choose the best combination of H curves from T . In many applications the number of curves of the raw template can be huge, thus we have developed a greedy algorithm which does not test all the possible combinations of curves but it is based on an iterative scheme. At each step, we divide T in subsets of Z curves, we look for the best one and we remove it from the raw template T saving it in a new template. This process is repeated until the new template has H curves, where H is the average number of curves in the population and Z is a submultiple of H. At each iteration, the best subset of Z curves is selected minimising:

arg min k ||T - (k+1)Z h=(kZ)+1 F h || 2 W * k = 0, ..., (H/Z) -1 (2.22)
after having randomized the order of the curves in T . The new template is the initial template used in the atlas construction.

Since the number of points per curve and the number of curves in the bundles can be dierent among the population, we have implemented this algorithm with the metric of varifolds, but any metric could be used.

An example of template initialisation applied to three brain ber bundles is shown in Fig. 2.1. The last row represents the initial templates used in the following experiments.

Multi-population atlas construction

In the previous subsections we assumed that all subjects belonged to the same group. In a clinical application, a subject can be characterised by a label which can be for instance his/her clinical status. Depending on this label, subjects can be 2.2. Methods divided into dierent groups. Here we propose to take advantage of this information computing a multi-population atlas construction.

Let L be the set of labels that can be given to a subject shape complex S i . We assume that L has only two elements for simplicity purpose: L = {c, p}, they can be interpreted as control and patient. The extension to more than two labels is straightforward but we will not consider it here. With this notation, we can rewrite Eq.2.1 for both groups as:

S c i = φ c i (T ) + ε c i S p i = φ p i (T ) + ε p i (2.23)
where we assume that the control points c and the template T are shared between the two populations. Instead, the initial momenta α 0i are divided into two classes which can be modelled as distinct Gaussian distributions:

α c i0 ∼ N (µ c , Γ c α ) α p i0 ∼ N (µ p , Γ p α ) (2.24)
We can quantitatively compare α c i0 with α p i0 since they are dened on the same set of control points and the groups share the same template. Moreover, in order to avoid identiability issues about the template, we impose that the sum of their Chapter 2. A Bayesian framework for multi-object atlas construction means µ c and µ p is equal to zero obtaining: µ c =µ and µ p =-µ.

We can then assume either that the two groups share the same variability (Γ c α =Γ p α =Γ α ) or that their within-group variations are dierent (Γ c α = Γ p α ). The computations for both cases can be found in 2.A.1. Using the second assumption, we could test the hypothesis of equality between the two covariance matrices or quantitatively compare their respective modes. On the contrary, if we ignored the information of the label computing a single-population atlas, we would implicitly assume that dierent groups belong to the same population and that they share the same variability. If we adopted another strategy by constructing an atlas for each group, we would obtain dierent templates and sets of control points. Thus, we could neither quantitatively compare the related populations of initial momenta nor directly compare the resulting templates.

The multi-population extension is computationally feasible only when using the proposed Bayesian framework. As previously shown, we use the log-likelihood of the initial momenta as regularity term of the deformations. Dividing the initial momenta into two dierent groups entails the addition of another regularity term into the cost function. This new term should be balanced by a trade-o parameter tuned using cross-validation. Increasing the number of groups would also augment the number of regularity terms and balancing values to tune. With the proposed statistical setting, there is no need to tune these weights since they are automatically taken into account in the estimates of the covariance matrices.

Experiments

In this section, we rst present the dataset used in the following experiments. Then, we evaluate the robustness of our Bayesian framework with respect to the values of the hyperparameters and to the number of points of the varifold grid. Furthermore, we compare these results with the ones of a previous atlas procedure without automatic estimates. Eventually, we present how the proposed method can be used to highlight morphological dierences between two groups using a multi-population atlas construction approach.

Materials

The dataset used throughout this chapter contains 40 subjects: 20 controls and 20 patients subject to Gilles de la Tourette syndrome. This neuropsychiatric disorder is thought to be associated with dysfunctions of the cortico-striato-pallido-thalamic circuits which are composed of sub-cortical structures linked to the cortical surface by ber bundles. We will investigate the morphological characteristics of these structures which distinguish the group of patients from the one of controls and which might be due to the Gilles de la Tourette syndrome.

For each subject, we consider three sub-cortical structures (left caudate, left putamen and left thalamus) and the ber bundles connecting them to the left hemi-

Robustness with respect to the hyperparameters

We evaluate here the robustness of the proposed algorithm with respect to the hyperparameter values {w j , P j , w α }. We do not consider P α since, as explained in Sec.2.2.4, we x it to P α = K -1 V . We compute 18 dierent atlases changing every time only one of the hyperparameters and keeping xed the others at a certain value.

Both w j and w α are previously normalised in order to use the same range of values, this means: w ε = w ε N Λ and w α = w α N giving as result:

Γα = 1 N N i=1 (α i0 )(α i0 ) T + w α P T α (1 + w α ) σ2 j = 1 N Λ j N i=1 ||S ij -φ i (T j )|| 2 W * + w j P j (1 + w j ) (2.25)
In order to understand a plausible range of testing values, we compute the resid-

uals 1 N Λ j N i=1 ||S ij -φ i (T j )|| 2
W * using the aforementioned dataset and initial tem- plates and we notice that the maximum value considering either only the sub-cortical structures or only the ber bundles is never above 10 and 100 respectively. This means that the product w j P j should not surpass this limit, otherwise a reduction of the data-term would impact very little the estimate of σ2 j and consequently the minimization of the cost function. Noting that, we choose as xed value 10 -2 and as testing range for all hyperparameters the interval [10 -2 -10 3 ] for the sub-cortical structures and [10 -2 -10 4 ] for the ber bundles.

In order to test the robustness of the results we compute the norm of the difference between the resulting template complex of all 18 atlases and a reference template complex. Moreover, we analyse also the norm of the residuals obtained at the end of the atlas constructions.

To show the importance of an automatic estimate of σ 2 j and Γ α , we compare our results with the ones obtained using a previous atlas procedure (Fixed) where both σ 2 j and Γ α are xed. In order to be coherent with the proposed Bayesian algorithm we set Γ α = K -1 V , which is the value of our prior P α and we let σ 2 j change along the same range of values of the hyperparameters. This means building 6 other atlases Chapter 2. A Bayesian framework for multi-object atlas construction changing only the value of σ 2 j . In all 24 atlas procedures, the other parameters, like λ W and λ V , are kept equal. Results using only left caudates or left caudate bundles of the 20 controls are shown in Fig. 2.2 where we have used as reference template the one obtained using the Fixed method with σ 2 j =10 -1 .

In Fig. 2.3 we compare the robustness of the results optimizing simultaneously two structures in each atlas construction: left caudate and left caudate bundle.

This means that we need to x two σ 2 j values with the Fixed method and 5 hyper- parameters in the Bayesian framework (w α , two w j and two P j ). Using as range of values 2.2.

The reference template is the one obtained at the end of the 9-th case using the Fixed method.

It is possible to conclude that the Fixed method leads to much more variable results than using Bayesian priors. The choice of the hyperparameters in our Bayesian framework is therefore easier than the one of σ 2 j . We decide to set them at: ω j =0.01Λ j N and P j =

0.05R 0 j ω j
where R 0 j is the initial data-term of structure j.

This brings to σ 2

j = R j +0.05R 0 j N Λ j (1+0.01)
which means that the minimum value of σ 2

j is equal to about 5% of R 0 j N Λ j
. This choice reduces the risk of over-tting since a value of R j smaller than

0.05R 0 j N Λ j
would almost not aect the estimate of σ 2 j . We also x ω α to 10 -2 . This makes negligible the contribution of P α to the estimate of Γα , especially when N is big, but it still permits to invert Γα .

Robustness with respect to the number of points Λ of the varifold grid

Another parameter of our algorithm is the number of points Λ of the varifold grid where template and shape complexes are projected to. We decided to make it dependent on λ W which is the bandwidth of the varifold kernel, as explained in Sec.2.2.2. In Fig. 2.4 we show that the results are not so sensitive to a change of the number of grid points and that our choice seems a good one. We compute 9 dierent atlases using the left-caudate surfaces of the 20 controls. All atlas constructions share the same set of parameters except for Λ which is tested in the range: [1-10 5 ].

As before we analyse the norm of the dierence between a reference template (Λ=1)

and the estimated templates at the end of the atlas procedures. We also show the norm of the residuals and the logarithm of the determinant of the covariance matrix of the initial momenta (log(|Γ α |)). As it is possible to notice there is a minimum value of the number of grid points ( 1000) above which the results are very similar

Experiments

and also more satisfactory.

In order to explain these results, it can be noticed that replacing σ 2 j with σ2 j into Eq.2.9 one obtains a new data-term which is equal to (w + ΛN ) log(σ 2 ). When Λ is small, the reduction of the residual, and consequently of log(σ 2 ), is not enough amplied by the multiplication with Λ to compensate a great increase in the regularity term. Thus, the estimated initial momenta keep a norm close to zero throughout the atlas construction. The deformations of the template are therefore insignicant and the covariance matrix Γ α depends only on the prior which explains why log(|Γ α |)

does not vary. The number of grid points obtained with the proposed method is equal to: 3456.

Chapter 2. A Bayesian framework for multi-object atlas construction

Experiments

Population dierences

In all the following atlas constructions we use the same parameters. All structures, except for the caudate, have a varifold parameter λ W equal to 2mm. For the caudate we choose a value of 3mm. About the bandwidth of the dieomorphic kernel K V , we use a value of 7mm with a consequent number of control points equal to 1080.

The maximum number of iterations is 120 and all the computations are performed on a Intel Xeon, 32 cores, CPU E5-2650, 2.60GHz with a graphic card NVIDIA Quadro 5000. All shape complexes are previously rigidly registered to a reference shape complex.

Multi-population Atlas

(Γ c α =Γ p α =Γ α )
We compute a multi-population atlas us- ing the two groups of controls and patients. The template complex and control points are shared between the two groups. Initial momenta are assumed to follow two Gaussian distributions with opposite mean and equal covariance matrix. The overall processing time is about 100 hours based on a CUDA implementation.

In Fig. 2.5 we show a scheme of the process. On the left there are the initial template complex and the initial set of control points which is initialised as an equally-distanced grid covering the entire ambient space. The distance is the bandwidth of the dieomorphic kernel: λ V . Momenta α i0 are initially all set to zero.

At the end of the atlas construction we obtain the nal template, the updated set of control points, the subject-specic sets of initial momenta, the average µ and the covariance matrix Γ α . The template shows the common shape features of both groups. At the top and at the bottom we present the average initial momenta for each group: µ and -µ. Then, on the right, we show the nal common template deformed using µ and -µ. The two complexes represent the anatomical congurations typical of each group and they can be directly compared since they both stem from the nal template complex. We take advantage of that by computing the absolute value of the dierence between the displacements from the nal template along the two average directions. This is shown in Fig. 2.6 for sub-cortical structures and in Fig. 2.7 for ber bundles. From the rst gure it is clear that the main differences are in the dorso-lateral part of the three sub-cortical structures, especially for the caudate. About the bundles, the dierences are mainly in the central part of the caudate-cortico and putamen-cortico bundles. Chapter 2. A Bayesian framework for multi-object atlas construction Multi-population Atlas (Γ c α = Γ p α ) Here we assume that the two groups of pa- tients and controls do not share the same covariance matrix. To test this assumption, we use the test proposed in Srivastava and Yanagihara (2010). The p-value (0.3461) is too big to reject the null hypothesis: Γ c α = Γ p α . This might be due to the fact that the sample size (20 subjects per group) is too small with respect to the number of deformation parameters (3240). Another possible explication could be that the two groups share most of their morphological variability which means that the rst modes of variation of both covariance matrices are similar. To test this hypothesis, we compute a Principal Component Analysis (PCA) for both covariance matrices and we calculate the angles between the modes of variation of the group of controls with the ones of the group of patients. The results for the rst 5 modes are shown in Table 2.1. As it is possible to notice, the rst mode of the group of controls is almost parallel to the rst one of the group of patients. This means that they produce similar morphological changes of the template complex. Moreover, they explain almost 35% and 45% respectively of the total variability of their groups, see Fig. 2.8. This might explain why the previous test fails in rejecting the hypothesis of equality between the two covariance matrices.

Experiments

The rst modes of both groups indicate the shared morphological variability between controls and patients. The second modes should instead produce shape variations which are more characteristic of each group (angle ∼ 60 • ) and which contribute also substantially (∼ 10%) to the total variability. In Fig. 2.9 we show the second mode of variation of both groups. Colors refer to the displacement from the nal common template shown in gray in the middle of the gure. It is possible to notice that the dorso-lateral part of the putamen varies almost in the same way between the two groups. Instead, the posterior-lateral part of the caudate and the ventral part of the thalamus present morphological variations which are dierent between the two groups. Chapter 2. A Bayesian framework for multi-object atlas construction

Experiments

Discussion and Conclusion

This work provides a Bayesian framework to embed a multi-object atlas construction into a statistical setting. It is general and it can be applied to any parametric deformation framework and to all shape models which permit to dene probability density functions. It allows to automatically estimate important balancing parameters which were xed by the user in previous methods, namely the noise variances of every structure. We have demonstrated that these parameters inuence much more the results than the hyper-parameters introduced with the proposed Bayesian priors. This statistical setting makes therefore multi-object atlas constructions more feasible and reproducible since the user is not obliged to tune or x balancing weights which grow with the number of analysed structures.

The proposed method allows also to estimate the covariance matrix of the deformation parameters. It is used for the statistical analysis of the morphological variability of the template complex within the population under study. In previ- A standard solution is to regularize it a posteriori using an identity matrix multiplied by a small scalar. This strategy is not satisfactory since it means that the estimate of the deformation parameters is based on a degenerate sample covariance matrix or on a wrong covariance matrix. Moreover, the regularization seems also too simplistic. In our approach, we directly estimate a well-conditioned covariance matrix using a more natural and coherent regularisation term given by the kernel of the RKHS to which belongs the vector eld used to compute the dieomorphisms.

Furthermore, the estimate of the covariance matrix does not change if the regularity term with its related prior terms are multiplied by a scalar. This means that we automatically take into account also the estimate of the so-called regularity level.

We have also extended the proposed statistical setting to a multi-population strategy where we take into consideration the aliation of a subject to a certain group dened, for instance, by its clinical diagnosis. This permits to employ more complicated but also more pertinent models at no additional cost, meaning without xing or tuning new parameters. One of the main advantages is to estimate dierent populations of initial momenta and not only one. This permits to quantitatively compare them since they are dened on the same set of control points and the groups share the same template. Moreover, this extension permits also to estimate a dierent covariance matrix for every group. We can therefore quantitatively compare them and their modes. It is interesting to notice that one could also compute the sample covariance matrices of dierent groups at the end of a single-population atlas and then compare their modes. But this would not be consistent with the fact that in a single-population atlas the initial momenta of dierent groups are estimated as belonging to the same Gaussian distribution.

We propose to model both curve and surface meshes as Gaussian random var-

Discussion and Conclusion

ifolds. They can be seen therefore as instances of the same mathematical object.

Since the space of varifolds is of innite dimension, we create a grid for every structure j on which we project both shapes and template. We decide to make the number of grid points Λ j dependent on the varifold bandwidth λ j W . This choice does not inuence the results in a single-object atlas as demonstrated in Sec.2.3.3.

In a multi-object analysis, it makes the contribution to the cost function (Eq.2.9) of structure j even more dependent on λ j W . In fact, the contribution of structure j depends mainly on the residual and on the value of Λ j which are both dened by λ j W . On the one hand, this choice simplies the atlas construction since the user has to x only one parameter, on the other hand it makes the choice of λ j W even more crucial. In future works, we will investigate how to automatically estimate λ j W . Another important parameter xed by the user is the dieomorphic bandwidth λ V . It is chosen based on the desired registration accuracy and it denes the number of control points. Its estimate, together with the one of λ j W , would make the atlas construction completely automatic whatever the number or kind of meshes. A possible strategy could be to select a range of suitable values of λ V and then use sparse multi-scale methods [START_REF] Sommer | Sparse Multi-Scale Dieomorphic Registration: The Kernel Bundle Framework[END_REF]. Otherwise, one could estimate the best deformation modules from a dictionary using sparse techniques in order to disentangle the single complicated dieomorphism into interpretable transformations [START_REF] Gris | A Sub-Riemannian Modular Approach for Dieomorphic Deformations[END_REF]. This would augment the computational load and execution time but it would also make the analysis more objective.

Another improvement for our method would be the employment of sampling algorithms like MCMC to ensure the convergence of the EM scheme whatever the quality of the data. As already said, these methods require a great computational burden. A possible solution would be to use a GPU implementation based on Multi-Graphics Processing Units (Multi-GPU) as in [START_REF] Ha | Fast Parallel Unbiased Dieomorphic Atlas Construction on Multi-graphics Processing Units[END_REF]. The authors showed that a Multi-GPU implementation gives a computational gain up to sixty times faster than a single CPU implementation.

We tested the proposed algorithm by comparing two populations of shape complexes, one of controls and one of Gilles de la Tourette patients, obtained from MR images of the brain and consisting of three sub-cortical structures and the ber bundles connecting them to the cortical surface. Results proved the eectiveness of our method in detecting morphological dierences between the two populations. They could indicate atypical connections resulting from abnormal brain development due to Gilles de la Tourette syndrome. Future works will aim to conrm this study by adding more structures as the cortical surface and possibly more subjects. Moreover, we plan to apply our algorithm to other human organs such as the heart or the liver.

Chapter 2. A Bayesian framework for multi-object atlas construction 2.A Appendix

2.A.1 Multi-population atlas construction

Here we derive the changes of the cost function (Eq.2.9) and relative gradients for the multi-population extension of the atlas construction.

2.A.1.1 Γ c α =Γ p α =Γ α
Let N c and N p be the number of subjects of the two groups and N =N c +N p . Moreover, we dene α c i0 ∼ N (µ, Γ α ) and α p i0 ∼ N (-µ, Γ α ). Thus, Eq.2.9 becomes:

M j=1 Nc i=1 1 2σ 2 j ||S c ij -φ c i (T j )|| 2 + P j w j N + 1 2 Nc i=1 (α c i0 -µ) T (Γ α ) -1 (α c i0 -µ) + M j=1 Np i=1 1 2σ 2 j ||S p ij -φ p i (T j )|| 2 + P j w j N + 1 2 Np i=1 (α p i0 + µ) T (Γ α ) -1 (α p i0 -µ) + M j=1 1 2 (w j + Λ j N ) log(σ 2 j ) + 1 2 (w α + N ) log(|Γ α |) + w α 2 tr((Γ α ) -1 P α ) (2.26)
where the norm ||•|| can be both the L 2 -norm and the varifold one. The gradients with respect to T , c 0 and {σ 2 j } are exactly the same as in Eq.2.17 and Eq.2.16. The one with respect to α i0 depends on the group and it is equal to:

∇ α c s0 E = M j=1 1 2σ 2 j ∇ α c s0 D c sj +Γ -1 α (α c s0 -µ) ∇ α p s0 E = M j=1 1 2σ 2 j ∇ α p s0 D p sj +Γ -1 α (α p s0 +µ) (2.27)
There is a closed form solution for both µ and Γ α which is equal to:

Γα = Nc i=1 (α c i0 -µ)(α c i0 -µ) T + Np i=1 (α p i0 + µ)(α p i0 + µ) T + w α P T α (N + w α ) μ = Nc i=1 α c i0 - Np i=1 α p i0 N (2.28) 2.A.1.2 Γ c α = Γ p α
As before, let N c and N p be the number of subjects of the two groups and N =N c +N p . Moreover, we dene α c i0 ∼ N (µ, Γ c α ) and α p i0 ∼ N (-µ, Γ p α ). The priors on Γ p α and

2.A. Appendix

Γ p α are both equal to W -1 (P α , w α ). Thus, Eq.2.9 becomes:

M j=1 Nc i=1 1 2σ 2 j ||S c ij -φ c i (T j )|| 2 + P j w j N + 1 2 Nc i=1 (α c i0 -µ) T (Γ c α ) -1 (α c i0 -µ) + M j=1 Np i=1 1 2σ 2 j ||S p ij -φ p i (T j )|| 2 + P j w j N + 1 2 Np i=1 (α p i0 + µ) T (Γ p α ) -1 (α p i0 -µ) + M j=1 1 2 (w j + Λ j N ) log(σ 2 j ) + 1 2 (w α + N c ) log(|Γ c α |) + w α 2 tr((Γ c α ) -1 P α )+ + 1 2 (w α + N p ) log(|Γ p α |) + w α 2 tr((Γ p α ) -1 P α ) (2.29)
where the norm ||•|| can be both the L 2 -norm and the varifold one. The gradients with respect to T , c 0 and {σ 2 j } are exactly the same as in Eq.2.17 and Eq.2.16. The one with respect α i0 depends on the covariance matrix of the group and its formulation is equal to the one in Eq.2.27. Furthermore, there is a closed form solution for both µ, Γ c α and Γ p α . The last two have the same expression of Eq.2.16 with a sum over N c and N p . Instead, the optimal value for the average momenta is :

μ = Nc i=1 (Γ c α ) -1 α c i0 - Np i=1 (Γ p α ) -1 α p i0 N c (Γ c α ) -1 + N p (Γ p α ) -1
(2.30)

2.A.2 Robustness Analysis

We present in Moreover, the great quantity of streamlines might impede understanding the pattern of a bundle making therefore dicult its interpretation.

In this chapter, we propose to approximate a ber bundle with a parsimonious representation of weighted streamline prototypes. We exploit the fact that many streamlines starting from seeds in the same voxel or in neighbour voxels share the same pathway and ending area. We approximate these streamlines with one of them, called prototype. We use a computational model for both streamlines and prototypes characterised by an explicit and easily computable metric. This allows us to control the approximation error and to select the streamlines which minimize it as prototypes.

Related Work

In the last years, there has been a great eort to compactly represent a ber bundle. A pragmatic strategy is to randomly choose a subset of the streamlines. The

Related Work

sampling is not driven by the minimisation of an approximation error and it is not possible to control the selection of the streamlines. This can cause the loss of the smallest fascicles of the bundle which might be important for the purpose of the study. More sophisticated solutions have been proposed and they can be separated into two categories. The rst group gathers the computational methods which compactly parametrise the single streamlines. The second category pools instead the strategies focused on concisely representing the entire ber bundle.

Compact representation of streamlines

Streamlines are composed of contiguous variable-length segments whose number might also vary among bers. Dierent computational models have been proposed

whose goal is to concisely parametrise a streamline. A widely employed method to easily compare and average streamlines but it is not a geometric primitive and it is therefore dicult to employ in multi-object registrations together with other geometric primitives such as 3D surface meshes.

Compact representation of the whole bundle

The second category is composed of methods which approximate the entire ber Moreover, the framework of currents does not take into account the extremities of the streamlines. This prevents the analysis of the structural connectivity, namely the areas of the gray matter connected by the bundle.

Our contribution

In this chapter, we propose to approximate any ber bundle with a set of weighted prototypes. Prototypes are chosen among the streamlines and they represent ensembles of similar bers. Their weights are related to the number of streamlines approximated. Both prototypes and streamlines are modelled as weighted currents, an extension of the framework of currents. This computational model takes into consideration both the pathway of the ber/prototype and the anatomical location of its extremities. Two bers/prototypes are therefore considered similar if their endpoints are close to each other and if their trajectories are similar. The space of weighted currents is a vector space with an explicit and easily computable metric.

This permits to approximate a bundle of streamlines controlling the approximation error. Moreover, the framework of weighted currents inherits from the one of usual The chapter is organised as follows. In Section 3.4 we rst present the framework of weighted currents and then the dierent steps of the proposed approximation scheme. After that, we demonstrate the eectiveness of our algorithm on deterministic and probabilistic ber bundles from both a qualitative and quantitative point of view. The last two sections are about the discussion of the results and conclusions.

Method

Weighted currents

The framework of weighted currents is an extension of the one of currents Vaillant and Glaunès (2005). A streamline is considered as a set of disconnected oriented points which are weighted by the spatial coordinates of the streamline extremities.

This means that every oriented point encodes not only the local orientation of a streamline but also its connectivity. It is an adaptation of the framework of functional currents Charon and Trouvé (2013a).

A streamline X is a polygonal line of N segments which is assumed to be an oriented and rectiable curve in R 3 . The coordinates of the two extremities f a and f b are two 3D vectors dened in the space Q=R 3 xR 3 . The ber X is modelled as a 1-weighted current C X via a line integral of a vector eld ω:

C X (ω) = X w(x, f a , f b ) T α(x)dx ≈ N i=1 ω T (x i ,f a ,f b ) (α i ) ≈ N i=1 δ α i (x i ,f a ,f b ) (ω) (3.1)
where x i and α i are respectively the center and the tangent vector of segment i which is approximated by a Dirac delta weighted current δ. We take advantage of the fact that the tangent vectors α(x) are constant over every segment and it can be shown that the approximation error tends to zero as the sampling becomes more accurate i.e. the length of the segments decreases [START_REF] Glaunès | Transport par diéomorphismes de points, de mesures et de courants pour la comparaison de formes et l'anatomie numérique[END_REF]. The vector eld ω belongs to a reproducing kernel Hilbert space (RKHS) W dened on the product space R 3 xQ. The space of weighted currents is a continuous linear form on W and every weighted current C X belongs to its dual space W * . A natural way to build a kernel K associated to the product space W is as tensor product of two kernels dened separately in R 3 (K g ) and in Q (K f ): K=K g ⊗ K f . Since even Q is a product space, K f is also dened as a tensor product between two kernels K a and K b . Thus, the kernel K results: K((x, f a , f b ), (y, t a , t b ))=K g (x, y)K a (f a , t a )K b (f b , t b ). All kernels K g , K a and K b are dened as Gaussian and they are parametrised by their bandwidths λ g , λ a and λ b . Using these kernels, the inner product in the frame-Chapter 3. Parsimonious approximation for white matter ber bundles based on weighted currents work of weighted currents between two Diracs is dened as:

δ α (x,f a ,f b ) , δ β (y,t a ,t b ) W * = K a (f a , t a ) K b (f b , t b ) (α T K g (x, y)β)
. By linearity, the inner product between two streamlines X and

Y (C Y (ω) ≈ M j=1 ω (y j ,t a ,t b ) (β j )) is: C X , C Y W * = K a (f a , t a )K b (f b , t b ) N i=1 M j=1 α T i K g (x i , y j )β j = exp -||f a -t a || 2 λ 2 a exp -||f b -t b || 2 λ 2 b (3.2) N i=1 M j=1 exp -||x i -y j || 2 λ 2 g α T i β j
The framed part would be the inner product between X and Y if modelled as usual currents. It measures overall dierences between the geometry of their trajectories. The two other terms take into account the dierences in the location of the streamline extremities. As shown in Fig. 3.1, the framework of currents is almost blind to a change of the positions of the end-points. Even if the extremities of the two streamlines are far from each other, with respect to the kernels' bandwidths, the angle (cos

-1 C X ,C Y W * |C X | W * |C Y | W *
) between the two streamlines varies by only few degrees. On the contrary, the streamlines become almost orthogonal in the framework of weighted currents when the extremities are apart. This makes the denition of similarity twofold in the framework of weighted currents. Two streamlines are considered similar if their pathways are alike, as in usual currents, but also if their endpoints are close to each other. In Fig. 3.2 we show the most similar bers to the red streamline in the framework of currents (green) and in the one of weighted currents (blue). The green bers share a similar pathway with the red one but some of them connect dierent anatomical areas. On the contrary, the blue bers are similar to the red streamline both in terms of geometry and connectivity.

As usual currents, the framework of weighted currents does not need point-topoint correspondence, except for the extremities which means that bers need to have a consistent orientation. Moreover, every streamline S i is considered as a vector in a Hilbert space. Thus, a ber bundle, which is the union of many bers B = ∪ i S i , is represented as a sum in this framework: C B = N i C S i . The dierence between two streamlines is dened as their sum with the orientation of the second ber inverted. If two bers are equal, their dierence cancels out. Furthermore, it is also possible to compute the average weighted current S of a ber bundle as:

C S = 1 N N i C S i .
Given the inner product dened in Eq.3.2, the squared norm of the dierence between two bundles in the framework of weighted currents. If the two streamlines X and Y are similar in this framework, their sum can be well approximated by C B =2C X or C B =2C Y . This is crucial for the scope of this chapter since it allows to represent an ensemble of streamlines with a weighted prototype where the weight is related to the number of streamlines approximated. In the previous example both X and Y could be chosen as prototype and the weight would be 2. A weighted prototype can be visualised as a tube where the streamline chosen as prototype is the central axis and the constant radius is proportional to the weight (see Fig. 3.3). In the following, we will describe how to use this idea to approximate a complex bundle stemming from a tractography algorithm. We will also assume that both streamlines and bundles are modelled as weighted currents writing simply S (resp. B) instead than C S (resp. C B ).

C B = N i=1 C S i and C B = M p=1 C S p is equal to: ||C B -C B || 2 W * = N i=1 N j=1 C S i , C S j W * + M p=1 M q=1 C S p , C S q W * -2 N i=1 M p=1 C S i , C S p W * . A bundle B composed of two streamlines X and Y is modelled as C B = C X +C Y 3.4. Method

Approximation scheme

The goal of the proposed approximation scheme is to represent a ber bundle B with a set of weighted prototypes {τ k P k }. It is based on a greedy approach where we rst subdivide the bundle into fascicles and then select the prototypes in each fascicle independently.

Fascicles detection A fascicle is a group of streamlines which are considered similar in the framework of weighted currents, namely they have a similar pathway Chapter 3. Parsimonious approximation for white matter ber bundles based on weighted currents

Q = N F F =1    i∈F S i 2 W * j / ∈F S j 2 W * -( i∈F j / ∈F S i , S j W * ) 2    (3.3)
where F is a fascicle, N F is the number of fascicles and it is constrained by

N F F =1 i∈F S i = B.
In the simple case of N F =2, Eq.3.3 can be rewritten as: Prototype Streamline Selection (PSS) Once dened the fascicles, a PSS is performed in each fascicle independently. We propose an iterative algorithm in the spirit of orthogonal matching pursuit [START_REF] Tropp | Signal Recovery From Random Measurements Via Orthogonal Matching Pursuit[END_REF]. Let F be a fascicle with L streamlines modelled as weighted currents, the rst prototype P 1 is chosen as the streamline minimising the residual squared error, namely:

Q = S1 2 W * S2 2 W * -S1 , S2
P 1 = arg min S i ||F -τ 1 S i || 2 W * .
Since the space of weighted currents is a vector space we can easily minimize this cost function. The optimal weight is:

τ 1 = F,P 1 W * ||P 1 || 2 W *
and the prototype is:

P 1 = arg max S i F, S i ||S i || W * 2 W * = arg max S i L 2 S, S i ||S i || W * 2 W * with
i=1,...,L. The prototype is therefore the most parallel streamline to the average S of the fascicle.

Once the rst prototype is selected, we remove from each streamline S i its orthogonal projection onto the prototype, resulting in the residual: r(

S i ) = S i -π(Si) = S i - S i ,P 1 W * P 1 ||P 1 || 2 W *
keeping therefore only the components of the streamlines orthogonal to the prototype P 1 . In this new space, we select the second prototype as:

P 2 = arg max r(S i ) r(F ), r(S i ) r(S i ) W * 2 W * .
We iterate this process until:

||F -K k=1 τ k P k || W * ≤ γ||F || W *
where ||F || W * is the norm of the fascicle, K is the number of prototypes and γ indicates the required approximation level. At each iteration t, the set of weights {τ k } k=1,...,t is computed as the orthogonal projection of all the streamlines of F to the space spanned by the selected set of prototypes Chapter 3. Parsimonious approximation for white matter ber bundles based on weighted currents

{P k } k=1,...,t .
It is important to notice that all these computations are based on the Gram matrix Γ of the fascicle F which has size [LxL]. Thus, instead of computing directly r(S i ), we simply update Γ as: r(S i ), r(S j ) W * = S i , S j W * -

S i ,P W * S j ,P W * ||P || 2 W * = Γ (i,j) - Γ (i,P ) Γ (j,P ) ||Γ (P,P ) || 2 2
. A sketch of the algorithm can be found in Algorithm 1 where Γ (i,j) indicates the value of the matrix Γ at row i and column j, Γ (K,L) refers to the submatrix of Γ containing the K rows of the prototypes and all the L columns, Γ (K,K) is the square submatrix with the rows and columns of the K prototypes and 1 (L,1) is a L-dimensional column vector of ones.

Algorithm 1 Prototype Fiber Selection Input: Fascicle F with L streamlines S i : F = L i=1 S i ; threshold γ ; pre-computed matrix Γ ; K=1 1:

P 1 ← arg max S i F, S i ||S i || W * 2 W * 2: τ 1 ← F,P 1 W * ||P 1 || 2 W * 3: for i = 1 to L do 4: for j = 1 to L do 5: S i , S j W * ← S i , S j W * - S i ,P 1 W * S j ,P 1 W * ||P 1 || 2 W * 6:
end for 7: end for 8: while ||F -K k=1 τ k P k || W * ≤ γ||F || W * do 9: K ← K + 1 10:

P K ← arg max i (Γ (i,L) 1 (L,1) ) 2 Γ (i,i) 11: {τ k } k=1,...,K ← Γ (K,L) 1 (L,1)
Γ (K,K)

12:

for i = 1 to L do 13:

for j = 1 to L do 14:

Γ (i,j) ← Γ (i,j) - Γ (i,P K ) Γ (j,P K ) ||Γ (P K ,P K ) || 2 2 15:
end for 16:

end for 17: end while Output: {τ k }, {P k } k=1,...,K After selecting the prototypes of each fascicle independently, they are all gathered into a single bundle of prototypes B P . The weights are then recomputed by the orthogonal projection of the whole bundle B to the entire set of prototypes B P in order to retrieve the correct values also for the weights of the prototypes close to the boundary between two dierent fascicles. Moreover, before the PSS algorithm, we perform also an outlier detection step in every fascicle. The streamlines characterised by an average angle with the other streamlines between 88 • and 90 • are considered as outliers and discarded from the analysis. This permits to reduce the number of analysed streamlines, and consequently also the amount of computations, Remark Performing the PSS in each fascicle independently permits to distribute the computations to dierent processors, reducing therefore the computational time.

Method

Moreover, it allows also to reduce the chance to select as prototype a streamline which might be considered an outlier. This is explained in Fig. 3.4 where every dot represents a streamline of a bundle modelled as weighted current. We oversimplify this space assuming it is simply R 2 . In this space, it is likely that a tractography bundle has a multi-modal distribution, where every mode is a fascicle. If we wanted to approximate the whole bundle with a single prototype, it would be the most parallel streamline to the average of the bundle. In Fig. 3.4 we would choose the streamline represented by the green cross. This ber is far from almost all the other streamlines and it could be considered as an outlier. Instead, if we apply the same selection process in each fascicle independently we would obtain the three prototypes highlighted in red. These streamlines are more representative than the green ber and they better approximate the bundle. Obviously, an actual fascicle does not lie in a 2D space and therefore we may need more than one prototype to explain its variability. Chapter 3. Parsimonious approximation for white matter ber bundles based on weighted currents

Experiments and Results

In this section, we rst describe the dataset used in the following experiments and some technical details about the implementation of the proposed algorithm. Then, we present the approximation of two probabilistic bundles and we show that their structural connectivity is similar to the one of the original bundles. Furthermore, we qualitatively evaluate the eect of the parameters on the approximation of a deterministic bundle. Eventually, we assess quantitatively the performance of our algorithm showing that the registration between two approximated ber bundles is denitely faster than using the original bundles for the same registration accuracy.

Materials

We test the proposed approximation scheme on 25 subjects.

Numerical aspects

The parameters needed to be xed by the user are the bandwidths of the three kernels of weighted currents λ g , λ a and λ b and the approximation level γ. In the following, λ g refers to the kernel of usual currents, λ a to the end-point on the nuclei and λ b to the end-point on the cortical surface.

All experiments shown in this chapter are computed on a Intel Xeon, 8 cores, CPU E5-1620, 3.60 GHz with a memory of 16Gb and a graphic card NVIDIA Quadro K4000. The computational times for the approximations of the ber bundles shown in

Weighted prototypes representation

We present in Fig. 3.5 the weighted prototype approximations of two probabilistic bundles using γ= 0.13, λ g =7mm, λ a =5mm and λ b =10mm. The proposed representation preserves the global shape of the bundle and it approximates thoroughly all the fascicles. We show also the probability densities of the endpoints (i.e. structural connectivity) on the cortical surface computed using both the streamlines of the original bundles and the weighted prototypes. The two densities are very similar from a qualitative point of view and the Kolmogorov-Smirnov test fails to show statistically signicant dierences between them at the 5% level. Probability densities are computed using Gaussian kernels, taking into account the weights of the prototypes for the proposed approximation. Furthermore, we use the same parameters to approximate all the 150 bundles of our data-set and in any case the density of the endpoints is statistically signicantly dierent from the one of the original bundle.

The average compression ratios, 100(1-K/N) where N is the number of streamlines of the original bundles and K is the number of prototypes, are shown in Table 3.5.3 and 3.5.3 for the deterministic and probabilistic bundles respectively. These results

show that our algorithm leads to a much more compact representation of the bundle while preserving the overall structural connectivity, which may be used for further analysis.

Experiments and Results

Figure 3.6: Weighted prototype representations (in red) at dierent approximation levels of a deterministic cortico-putamen bundle (in blue). Each row is based on a dierent set of kernel bandwidths. The letters N and K refer respectively to the number of streamlines of the bundle and to the number of prototypes. The compression ratios are indicated in brackets.

Qualitative evaluation of the parameters inuence

In Fig. 3.6 we evaluate the inuence of the parameters of our algorithm on a deterministic cortico-putamen bundle. In the rst row we employ λ g =5mm, λ a =4mm, λ b =6mm and in the second row λ g =7mm, λ a =5mm, λ b =10mm. Every column corresponds to a dierent approximation level. It can be noticed that at γ=0.05 and γ=0.13 all fascicles are well approximated, whatever the set of parameters. At γ=0.25 and higher values of γ (not shown here) only the denser parts of the fascicles are kept. Moreover, the results based on the rst row use almost twice the prototypes than in the second row. This is expected since the values of the kernels are smaller and therefore the denition of similarity between two streamlines is more Chapter 3. Parsimonious approximation for white matter ber bundles based on weighted currents stringent. Thus, for a given γ, one needs more prototypes to approximate the same number of streamlines.

Furthermore, this inuences also the number of fascicles, as it can be seen from

Registration-based evaluation of the algorithm

Here we evaluate the impact of the proposed approximation scheme on the quality of a registration between two deterministic cortico-putamen bundles of dierent subjects. We use the dieomorphic transformation implemented in the software Deformetrica (www.deformetrica.org) with a kernel bandwidth equal to 10mm. The source bundle B S is the one shown in Fig. 3.6. First, we approximate both B S and the target bundle B T , composed of 25916 streamlines, at dierent approximation levels using λ g =7mm, λ a =5mm, λ b =10mm. Then, for each level, we register the approximation of B S onto the one of B T . We apply then the obtained deformation to the original bundle B S and we measure the residual error between the transformed original source bundle φ(B S ) and the original target bundle B T :

||φ(B S ) -B T || 2
W * in the framework of weighted currents. Ideally, we would compare the obtained residual errors with the one of the registration between the original ber bundles.

Unfortunately, the computational time would be too long (see Table 3.5.5) and therefore we decide to use the registrations between smaller sub-samples of the original ber bundles for comparison. In Fig. 3.8 we show the deformations of the original ber bundles using the deformations obtained with the approximation at γ=0.13 and with the sub-sample of 5000 streamlines. It is possible to notice that the results look very similar. This is conrmed in Table 3.5.5 where we show that the dierence between their residual errors is very small. We can conclude that the registration based on the approximation at γ=0.13 is as accurate as using the sub-sample of 5000 streamlines but 93 times faster! We present also the results for other approximation levels. With respect to γ=0.13, the registrations are either less accurate or slower and with a similar accuracy.

Experiments and Results

Registration

Chapter 3. Parsimonious approximation for white matter ber bundles based on weighted currents

Discussion and Conclusions

We presented an algorithm to approximate a ber bundle with a small and scattered set of weighted streamlines prototypes. We tested it on 150 bundles resulting from both deterministic and probabilistic tractography algorithms. The number of prototypes was on average 2% of the number of streamlines of the original bundles.

We showed that such a parsimonious representation preserves both the shape and the structural connectivity of the bundles.

The streamlines considered in this chapter have been truncated at the intersection between gray and white matter which is an area usually characterised by a low Signal to Noise Ratio. This means that the estimates of the end-points of the streamlines are not very robust. To account for this uncertainty, we use Gaussian kernels for measuring the dissimilarity between two streamlines. Theoretically, the greater the uncertainty, the greater the bandwidths of the kernels. In the proposed algorithm, these bandwidths (i.e.

(1 -γ)||F || W * ≤ || K k=1 τ k P k || W * ≤ (1 + γ)||F || W * .
This means that, using γ=0.13, the norm of the prototypes will be at least 87% of the norm of the fascicle at the end of the algorithm. Furthermore, we noticed that between γ=0.5 and γ=0.01 the number of prototypes grows exponentially. For instance, using the set of parameters of Fig. 3.5, the number of prototypes is 70 at γ=0.5, 367 at γ=0.13 and 2307 at γ=0.01. This means that we need less prototypes to reduce γ from 0.5 to 0.13 than from 0.13 to 0.01. This is because the rst prototypes approximate the parts of the fascicle with an higher density of streamlines (i.e. greater redundancy). Thus, their weights have a great value and few prototypes can explain a considerable percentage of the norm of the fascicle. Instead, between γ=0.1 and γ=0.01, most of the streamlines have already been approximated and every new prototype can explain only a few of the left bers. We found that a value of γ=0.13 results in a parsimonious representation which exhaustively approximates all the fascicles. An interesting improvement would be to automatically estimate the value of γ taking into consideration both the approximation error and the number of prototypes.

We showed also that our representation can approximate not only the central and more dense mass of a bundle, as usual currents, but also its smaller fascicles and extremities. We demonstrated its usefulness in the registration between two

3.A. Appendix

approximated bundles where we correctly match the whole pathway of the weighted prototypes, including their extremities (see Fig. 3.8). Furthermore, bundle registrations based on our approximation scheme present a lower computational time and memory footprint than using the original ber bundles. This makes thus possible population studies, like the atlas construction, based on multiple template-subject registrations, which would be infeasible with the original ber bundles.

3.A Appendix

3.A.1 Interpretation of τ

We show here that the weight τ of a prototype P is related to the number of bers approximated by P . Given a bundle B composed of 3 bers: B = 3 i=1 S i we want to approximate it with one prototype P . Let assume that S 1 is the prototype P (the reasoning does not change modifying the prototype), the value of its weight τ is 1 +

||S 2 || 2 W * ||P || 2 W * cos(P, S 2 ) W * + ||S 3 || 2 W * P 2 W *
cos(P, S 3 ) W * . This means that if S 1 is parallel to the other bers and their norms are similar, the value of τ will be about 3. Instead, if either S 2 or S 3 is orthogonal to S 1 , the prototype will not approximate that ber and the weight will be smaller than 3. This shows that τ is related to the number of bers approximated by the prototype or, more precisely, to how much similar is the prototype with respect to the other bers of the bundle. When dealing with more prototypes, every τ depends also on the inner product between the prototypes.

Since we project at each iteration all the streamlines onto the orthogonal space of the last estimated prototype, the inner product between prototypes should be small.

3.A.2 Modularity based on weighted currents

The denition of modularity in Blondel et al. ( 2008) is:

Q = N C c=1 W c m - S c 2m 2 (3.4)
where N C is the number of modules (fascicles In order to analyse the relative position between dierent structures, one needs to study the neural circuits as a whole. This requires a single framework where both grey and white matter components are studied together and at the same time.

), W c =1/2 i∈c S i 2 W * is the
Grey matter structures are usually analysed with T1-w scans whereas white matter tracts using dMRI scans. As already explained in the Introduction, we integrate the information from these two modalities by shifting the morphological analysis from images to complexes of 3D meshes, called shape complexes, representing the dierent components of the neural circuits. Structures of the gray matter such as cortical surface and basal ganglia are represented as closed surfaces segmented from T1-w scans whereas the neural projections of the white matter are modelled as bundles of 3D streamlines, called ber bundles, resulting from tractography algorithms applied on diusion scans.

Dierent strategies exist for the statistical shape analysis of a population of shape complexes, representing each one the structures of a subject. One which naturally allows the combination of dierent mesh types is based on the Grenander's pattern theory [START_REF] Grenandner | General pattern theory -A mathematical study of regular structures[END_REF]. Every shape complex is modelled as a deformation of a reference complex called template complex. The deformations are dened on the whole ambient space where the template complex lies and they put into correspondence the template complex with the homologous structures of the subjects. First examples used a xed template chosen a priori like the Talairach atlas [START_REF] Davatzikos | A computerized approach for morphological analysis of the corpus callosum[END_REF] or a randomly-selected subject [START_REF] Csernansky | Hippocampal morphometry in schizophrenia by high dimensional brain mapping[END_REF]. However, this kind of templates might have imperfections due to segmentation errors and they might not be centred with respect to the population under study, biasing consequently the whole analysis. Later works have thus proposed to automatically estimate the template [START_REF] Joshi | Landmark matching via large deformation dieomorphisms[END_REF][START_REF] Durrleman | Morphometry of anatomical shape complexes with dense deformations and sparse parameters[END_REF][START_REF] Qiu | Atlas Generation for Subcortical and Ventricular Structures With Its Applications in Shape Analysis[END_REF]. The joint estimate of the template complex and the deformations has been called atlas construction. The template complex captures the common morphological characteristics of the population. The set of deformations explain the variations in shape across the population.

Introduction

Deformations are dened as dieomorphisms of the whole ambient space which are smooth invertible transformations with smooth inverse. This kind of deformation preserves the anatomical organisation of the components of the template complex, namely they can not intersect, fold or shear. Moreover, since the transformations are dened locally and they can vary across dierent areas in the ambient space, it is also possible to capture the variations in relative position between separated structures. However, using a single dieomorphism we implicitly assume that the relative position between structures in contact with each other or close to each other does not change across subjects. This implies that a particular ber bundle should link the same areas of the cortical surface and basal ganglia in each subject within the population. This assumption precludes the study of changes in structural connectivity which could be caused by an abnormal brain development. In Fig. 4.1 we present a toy example composed of a template complex and a subject shape complex characterised by a dierent structural connectivity. A single dieomorphism could not put into correspondence all the structures and capture the dierences in structural connectivity.

Structural connectivity analysis is usually tackled by dividing the cortical surface and the sub-cortical nuclei in parcels based on anatomical, histological or functional schemes reproducible across subjects [START_REF] Craddock | Imaging human connectomes at the macroscale[END_REF]. Every parcel is then considered as a node of a graph and the number of streamlines connecting two nodes represents their weighted edge. Variability in structural connectivity across subjects can be dened by simply analysing independently the variations in each parcel or by using indexes and methods from the complex network theory (Bullmore and Sporns, 2009; [START_REF] Worbe | Altered structural connectivity of cortico-striatopallido-thalamic networks in Gilles de la Tourette syndrome[END_REF]. In both cases, the analysis depends on the chosen parcelling scheme and it does not take into consideration the morphological variability of the gray matter structures and the trajectories of the streamlines.

Here, we propose to join shape and structural connectivity analysis in a unied framework based on a double dieomorphic atlas construction. The template complex is warped towards every shape complex of the population using a cascade of two dieomorphisms. The rst dieomorphism acts only on the white matter of the template complex keeping xed the gray matter. It can be seen as a relative change of coordinates within the gray matter which is considered as a xed reference frame.

During this deformation the ber bundles are detached from the xed gray matter structures and re-attached in another location, capturing the variation in structural connectivity. These variations can be compared among subjects since they are all computed with respect to the same reference frame, namely the xed gray matter of the template complex. The second dieomorphism acts on the whole template complex, namely on both the resulting white matter and on the gray matter. This is a global change of coordinates which has the role to put into correspondence both gray and white matter structures of the template complex with the homologous ones of the subject's shape complex. The rst dieomorphism re-positions the white matter tracts within the gray matter such that all the components of the template complex can be correctly put into correspondence by the second diffeomorphism. The two dieomorphisms are optimised together minimising a single This is achieved by ensuring that the components of the displacement eld orthogonal to the motion boundaries vanish at the boundaries. In this way, every region deforms smoothly and independently from the others and discontinuities are allowed only at the boundaries. In our case, we are interested in studying the relative variation of one region, white matter, with respect to another one, gray matter. The aforementioned sliding registrations, if applied to the template complex shown in Fig. 4.1, would result in two independent deformations, one for the white matter and one for the grey matter, and it would be impossible to understand whether the changes of the white matter depend on a dierence in gray matter or in structural connectivity. Moreover, due to the sliding conditions, ber bundles should always 4.1. Introduction stay in contact with the grey matter surface during deformation. This would require a lot of deformation energy to move the ber bundle from the left gyrus to the right one and it would produce a very complicated deformation. Instead, with the proposed method the ber bundle is not constrained to slide on the cortex and therefore the resulting deformation would be simpler. Note that, we enforce that the template ber bundle will be near the grey matter at the end of the double dieomorphic deformation by modelling streamlines as weighted currents. As explained in Chapter 3, this is an extension of the framework of currents and two streamlines are considered similar not only if their pathways are alike but also if their endpoints are close to each other. This makes possible to match correctly also the extremities of two ber bundles and not only their central part as in usual currents. The template ber bundle is thus deformed in such a way that at the end of the second dieomorphism its extremities match the ones of the subject bundle and therefore they are also close to the grey matter structures.

In the proposed method, both dieomorphisms are parametrized using control points as presented in [START_REF] Durrleman | Optimal Data-Driven Sparse Parameterization of Dieomorphisms for Population Analysis[END_REF]. The number of control points is xed by the user and their positions are automatically adjusted close to the most variable parts of the template complex during the atlas construction.

In order to deal with the considerable amount of streamlines resulting from tractography algorithms, we rely on the approximation scheme introduced in Chapter This chapter is organised as follows. In Section 4.2, we rst present a double dieomorphic generative model which is integrated into a Bayesia framework for atlas construction. Afterwards, we recapitulate how to model surface meshes as varifolds (Charon and Trouvé, 2013b) and streamlines as weighted currents (Chapter 3) and how to dene probability densities in these innite-dimensional spaces.

Then, we present how to compute double dieomorphic deformations based on a control-points formulation of the LDDMM framework and how to integrate them in the proposed statistical setting. Furthermore, we use this method to compare the organisation and morphology of the cortico-putamen neural circuit of a population of healthy controls with the ones of GTS patients. Eventually, we show that the proposed double dieomorphic scheme discriminates better between controls and GTS patients than using a single dieomorphism.

Methods

Double Dieomorphic Generative Model

The proposed atlas construction is based on a generative statistical model. We assume that the population under study is composed of N subjects. Every subject i is characterised of M structures segmented from both structural and diusion images. We dene the mesh of structure j belonging to subject i as S ij . Every subject shape complex S i , dened for the moment in a generic way as the ensemble of all the meshes S ij of subject i, is modelled as a double deformation of a common template complex T plus a residual noise ε i . Both T and ε i are also dened as the ensembles of the templates T j and residuals ε ij . The rst deformation φ W acts only on the white matter of the template complex: T W . The second deformation φ All deforms both the resulting white matter φ W (T W ) and the gray matter of the template complex T G . This formulation derives from the forward model where we assume that all elements belong to an algebraic structure where addition is dened [START_REF] Allassonnière | Towards a Coherent Statistical Framework for Dense Deformable Template Estimation[END_REF][START_REF] Ma | Bayesian template estimation in computational anatomy[END_REF][START_REF] Durrleman | A forward model to build unbiased atlases from curves and surfaces[END_REF][START_REF] Durrleman | A statistical model of white matter ber bundles based on currents[END_REF] and it results:

S i = φ All i φ W i (T W ) ∪ T G + ε i (4.1)
The two deformations φ W i and φ All i , proper to subject i, are two dieomorphisms of the entire ambient space. They follow one another creating a cascade of dieomorphisms. The rst one φ W i , called also white dieomorphism, acts only on the white matter of the template complex. It represents a relative change of coordinates with respect to the gray matter of the template complex which is kept xed.

The second (or global) dieomorphism φ All i transforms the whole template complex, both white and gray matter, and it is a global change of coordinates which brings the template complex to the subject space.

Methods

The rst deformation φ W i makes the ber bundles slide onto the gray matter structures. The streamlines are detached and re-attached in another area such that the whole template complex, both gray and white matter, can be well matched to the subject shape complex by the second dieomorphism φ All i . The rst dieomorphism φ W i captures the dierence in structural connectivity, namely the variations in the connected areas of the gray matter. The second dieomorphism φ All i describes the global morphological changes of both white and gray matter.

It is important to notice that it is fundamental to rst deform the white matter of the template complex and then the gray matter in order to retrieve the main variations in structural connectivity within the population. In fact, the transformations due to the rst dieomorphisms {φ W i } are comparable across subjects since they are all computed with respect to the same reference frame, namely the xed gray matter of the template complex. If we changed the order, deforming rst the gray matter and then the white one, we would not be able to compare the variations in structural connectivity since the reference frame given by the gray matter would be specic to each subject. Furthermore, a question that naturally arises using the proposed method is about the uniqueness of the decomposition into two dieomorphisms in the regions containing only white matter structures. In these areas, the white matter could be deformed in two dierent but equivalent ways. A possible solution to obtain a unique decomposition is to choose a scale of deformation so that white matter objects are deformed by the second dieomorphism in a correlated way with respect to at least one gray matter structure.

Bayesian Atlas Construction

The goal of the atlas construction is to estimate the template complex T = T W ∪T G , the variations in structural connectivity within the population described by the ensemble of white dieomorphisms {φ W i } and the global morphological variations captured by the set of global dieomorphisms {φ All i }. Both dieomorphisms are parametrised by a set of parameters called respectively α W i and α All i specic to each subject i. We assume that these parameters follow a Gaussian distribution with zero mean and covariance matrix equal to Γ W α and Γ All α respectively:

α W i ∼ N (0, Γ W α ) p(α W i |Γ W α ) ∝ 1 |Γ W α | 1/2 exp - 1 2 (α W i ) T (Γ W α ) -1 α W i (4.2) α All i ∼ N (0, Γ All α ) p(α All i |Γ All α ) ∝ 1 |Γ All α | 1/2 exp - 1 2 (α All i ) T (Γ All α ) -1 α All i
Both distributions are completely described by their covariance matrix. Moreover, as usual in statistical learning, we assume that the residuals ε i follow a Gaussian distribution centred at 0. For now, we model all structures of the shape complexes S i and of the template complex T with landmarks. Every structure j has a Chapter 4. A deformation framework to unify morphometry and structural connectivity analysis number of landmarks equal to Λ j and the norm of the dierence between two meshes of the same structure is dened as the square root of the sum of squared dierences between pair of landmarks (L 2 -norm, || • || 2 ). Thus, the likelihoods of the residuals of white and gray matter structures modelled as landmarks are respectively:

ε W ij ∼ N (0, σ 2 j Id Λ j ) p(ε ij |σ 2 j ) ∝ 1 |σ 2 j | Λ j /2 exp - 1 2σ 2 j ||S ij -φ All i φ W i (T W j ) || 2 2 (4.3) ε G ij ∼ N (0, σ 2 j Id Λ j ) p(ε ij |σ 2 j ) ∝ 1 |σ 2 j | Λ j /2 exp - 1 2σ 2 j ||S ij -φ All i T G j || 2 2
In Sec.4.2.3 we will make clear how to adapt these equations when a mesh is modelled as varifold or weighted current. Whatever the model employed, the variance depends exclusively on the structure-dependent parameter σ 2 j which is also assumed to be a random variable. Moreover, from Eq.4.1 and Eq.4.3 it follows that all shapes {S ij } follow a Gaussian distribution:

S W ij ∼ N (φ All i φ W i (T W j ) , σ 2 j Id Λ j ) S G ij ∼ N (φ All i T G j , σ 2 j Id Λ j ) (4.4)
We can thus reformulate the goal of the atlas construction as estimating the template complex T , the two covariance matrices of the deformation parameters Γ W α and Γ All α and σ 2 j , knowing the shape complexes of the population {S ij } and assuming they follow a Gaussian distribution as in Eq.4.4. This can be achieved by maximizing the joint posterior distribution of T , σ 2 j , Γ W α and Γ All α , which are also modelled as random variables. Assuming independence between all random variables and considering {α W i } i=1,...,N and {α All i } i=1,...,N as hidden variables, it results:

{T * , Γ W * α , Γ All * α , σ 2 * j } = arg max T ,Γ W α ,Γ All α ,σ 2 j N i M j p(T j , Γ W α , Γ All α , σ 2 j , α W i , α All i , S ij )dα All i dα W i (4.5)
Not using priors for σ 2 j , Γ W α and Γ All α can produce degenerate estimates. Instead, as demonstrated in [START_REF] Allassonnière | Towards a Coherent Statistical Framework for Dense Deformable Template Estimation[END_REF], the introduction of inverse Wishart distributions makes possible to obtain good estimates even when the number of subjects is small. Their probability density functions are: 4.2. Methods σ 2 j ∼ W -1 (P j , w j ) p(σ 2 j ; P j , w j ) ∝ (σ 2 j ) -w j 2 exp -1 2

w j P j σ 2 j (4.6) Γ W α ∼ W -1 (P W α , w W α ) p(Γ W α ; P W α , w W α ) ∝ |Γ W α | -w W α 2 exp - 1 2 w W α T r((P W α ) T (Γ W α ) -1 ) (4.7) Γ All α ∼ W -1 (P All α , w All α ) p(Γ All α ; P All α , w All α ) ∝ |Γ All α | -w All α 2 exp - 1 2 w All α T r((P All α ) T (Γ All α ) -1 ) (4.8)
The scalars w j , P j , w W α and w All α are strictly positive and P W α and P All α are positive symmetric matrices. Since the maximization of Eq.4.5 is not tractable analytically, we use the EM (Expectation Maximization) algorithm where we approximate the conditional distribution of the E step with a Dirac distribution at its mode. Assuming that the template T has a non-informative prior distribution, it results:

M j=1 N i=1 1 2σ 2 j ||S ij -φ All i φ W i (T W j ) | 2 2 + P j w j N + M j=1 N i=1 1 2σ 2 j ||S ij -φ All i T G j || 2 2 + P j w j N + M j=1 1 2 (w j + Λ j N ) log(σ 2 j )+ 1 2 N i=1 (α All i ) T (Γ All α ) -1 α All i + (w All α + N ) 2 log(|Γ All α |) + w All α 2 tr((Γ All α ) -1 P All α )+ 1 2 N i=1 (α W i ) T (Γ W α ) -1 α W i + (w W α + N ) 2 log(|Γ W α |) + w W α 2 tr((Γ W α ) -1 P W α ) (4.9)
This equation represents the cost function of our algorithm. We use a gradient descent scheme to minimize it. The framed terms refer respectively to the dataterms and to the regularity terms of both dieomorphisms. The other terms are due to the use of the inverse Wishart prior distributions. The proposed statistical framework is general since it can be employed with any shape model, provided it is possible to dene probability density functions, and any parametric dieomorphic model.

Chapter 4. A deformation framework to unify morphometry and structural connectivity analysis

Shape models

The framework of landmarks has been extensively used to model shapes since the metric and its gradients are easy to compute and it is also simple to dene random variables. However, the denition of correspondences between 3D meshes is a tedious and dicult task which might also be impossible when dealing with ber bundles.

To this end, it is better to opt for correspondence-free shape models such as the framework of varifolds or currents. We propose to model gray matter structures as varifolds and white matter streamlines as weighted currents. Moreover, since ber bundles resulting from tractography algorithms might have a number of streamlines so big to make infeasible the atlas construction, we propose to approximate them with the parsimonious representation of weighted streamlines prototypes presented in Chapter 3.

In the following, we will rst recapitulate the framework of varifolds and weighted currents, then we will explain how to dene random variables based on these two shape models.

Varifolds

The framework of varifolds is the non-oriented extension of the one of currents [START_REF] Vaillant | Surface Matching via Currents[END_REF]. It can be used to model both surfaces and curves as in Chapter 2 but here we will use it only for surfaces. The denition of varifold for a 3D rectiable surface X is :

V X (ω) = X ω(p, ← -→ n(p))|n(p)| 2 dp ≈ L l=1 ω(p l , ← → n l )|n l | 2 (4.10)
where n(p) ∈ R 3 is the normal of X at the point p ∈ R 3 . Since the surfaces used throughout this chapter are meshes composed of a nite number of faces, it seems reasonable to approximate the integral in Eq.4.10 as a sum over the centres {p l } of the L faces of the surface as in Chapter 2. However, it is still an open problem the demonstration of the convergence of this approximation when the area of the faces decreases, i.e. the sampling becomes more accurate [START_REF] Charlier | The fshape framework for the variability analysis of functional shapes[END_REF]. The test eld ω ∈ W is function of both the position p and the unoriented unit normal vector ← -→ n(p). The set of unoriented unit vectors ← → n ∈ ← → S can be formally dened as the quotient of the unit sphere in R 3 by the two elements group {±Id 3 }. This means that two unoriented unit vectors which are superimposed (and thus parallel) but with a dierent orientation are considered equivalent. Eq.4.10 shows that the space of varifolds W * is a linear functional from W to R. This makes W * a vector space and therefore the union of meshes is equal to a sum in the space of varifolds.

But, contrary to currents, inverting the orientation of a mesh does not change the varifold representation.

The test space W is dened as a product space between the ambient space R 3 and the space of unoriented unit normal vectors

← → S , W : (R 3 x ← → S ). More

Methods

precisely, as proposed in (Charon and Trouvé, 2013b), one chooses W to be a Reproducible Kernel Hilbert Space (RKHS) whose kernel K W is dened as the tensor product between two kernels k g , Gaussian (continuous red line), and k n , Cauchy-Binet (dashed blue line), dened on R 3 and on ← → S respectively:

K W ((p, ← -→ n(p)), (q, ← -→ u(q))) = k g (p, q)k n ( ← -→ n(p), ← -→ u(q)
). Thanks to the reproducing property of the RKHS, the inner product between X and another surface Y modelled as varifold

(V Y ≈ H h=1 ω(q h , ← → u h )|u h | 2 ) results: V X , V Y W * = L l=1 H h=1 exp -||p l -q h || 2 2 λ 2 W n T l u h |n l | 2 |u h | 2 2 |n l | 2 |u h | 2 (4.11)
The only parameter is the standard deviation λ W of the Gaussian kernel. The distance between V X and V Y is:

||V X -V Y || 2 W * = V X , V X W * + V Y , V Y W * - 2 V X , V Y W * .
An important characteristic of this metric is the absence of correspondences: every normal of the rst surface is compared with all the normals of the other surface (Cauchy-Binet kernel) and their contribution is weighted by the distance between the centres (Gaussian kernel). Moreover, it is interesting to notice that the Cauchy-Binet kernel is equal to the squared cosine of the angle between the normals (i.e. n l , u

h 2 =n T l u h =|n l | 2 |u h | 2 cos(•))
. This makes the metric invariant to a change of orientation of some normals of the surfaces. Moreover, this property is particularly important in the estimation of the template since it permits to prevent the so called currents cancelling eect. This happens in surfaces modelled as currents when two faces with opposite orientation, like in a thin protrusion, cancel reciprocally their eect on the estimate of the template. For a more detailed discussion about varifolds the user is referred to Charon and Trouvé (2013b).

Weighted currents

The framework of weighted currents, presented in Chapter 3, is also an extension of the one of currents. It has been conceived for curves and particularly to model the streamlines resulting from tractography algorithms. An oriented 3D polygonal curve A composed of G segments is modelled as a 1-weighted current C A via a line integral of a square-integrable vector eld ω:

C A (ω) = A w(x, f a , f b ) T α(x)dx ≈ G g=1 w(x g , f a , f b ) T α g (4.12)
where x g ∈ R 3 and α g ∈ R 3 are respectively the centre and the tangent vector of segment g and, as previously for varifolds, we approximate the integral with a sum over the segments of A. The two 3D vectors f a and f b are the coordinates of the end-points of the curve and they are dened in the space Ψ = R 3 xR 3 . Similarly to varifolds, the vector eld ω belongs to a RKHS Q dened on the product space Chapter 4. A deformation framework to unify morphometry and structural connectivity analysis R 3 xΨ and the space of weighted currents Q * is a continuous linear form on Q. Thus, even in this case, every curve is considered as a vector in a Hilbert space and the union of several curves, such as a white matter ber bundle, is represented as a sum.

The kernel K Q associated to the RKHS Q is dened as:

K Q ((x, f a , f b ), (y, t a , t b )) = k g (x, y)k a (f a , t a )k b (f b , t b )
where the three kernels k g , k a , k b are Gaussian and dened in R 3 . As explained in Chapter 3, the inner product between A and another curve B modelled as weighted current (C B (ω)

≈ F f =1 w(y f , t a , t b ) T β f ) results in: C A , C B Q * = G g=1 F f =1 exp -||x g -y f || 2 2 λ 2 g α T g β f (4.13) exp -||f a -t a || 2 2 λ 2 a exp -||f b -t b || 2 2 λ 2 b
The continuous red line highlights the inner product between A and B if modelled as usual currents. It measures overall geometrical dierences between their pathways. The two other terms, underlined with dashed lines, constrain the denition of similarity. Two curves are considered similar if their pathways are alike, as usual currents, but also if their endpoints are close to each other. The inner product is parametrised by the three bandwidths: λ g , λ a , λ b . The distance between two curves is dened as:

||C A -C B || 2 Q * = C A -C B , C A -C B Q * .
As usual currents, this framework does not need point-to-point correspondences, except for the endpoints. This means that curves need to have a consistent orientation. This can be complicated when comparing whole-brain tractograms but here we need to compare only ber bundles extracted using two specic ROIs (Region Of Interest) dening the initial and last point respectively. This permits to easily ensure a consistent orientation between the streamlines of the same bundle and also among bundles of dierent subjects. Furthermore, the use of Gaussian kernels on the extremities of the streamlines permits to be more robust towards tractography errors at the boundary between grey and white matter which is an area usually characterised by a low SNR (Signal to Noise Ratio).

Varifolds and Currents random variables

In order to integrate the models of varifolds and weighted currents to the statistical framework previously presented, we need to dene random Gaussian variables in these two spaces. In both cases, they are dened as a linear map between every vector eld ω ∈ W (or ω ∈ Q) to a real random Gaussian variable G(ω) such that E[G(ω)] = 0 and, given two vector elds ω 1 and ω

2 , E[G(ω 1 ), G(ω 2 )] = ω 1 , ω 2 W (orQ)
. This means that the kernels K W and K Q dene the covariance structures of the Gaussian varifolds and weighted currents respectively. Since both of them are innite-dimensional, we dene for each structure j a nite dimensional space on which we project both the template T j and the shapes {S ij } i=1...N mod-4.2. Methods elled as varifolds or weighted currents. We follow the line of reasoning presented in [START_REF] Bibliography Durrleman | Statistical models of currents for measuring the variability of anatomical curves, surfaces and their evolution[END_REF] for the framework of currents.

Varifolds The nite dimensional space W * Λ j for varifolds is dened as the span of a set of delta Dirac varifolds: Span{δ (xu, ← → β k ) } where the points {x u } and the unoriented unit vectors { ← → β k } are constrained to belong to two predened grids, respectively Υ j

x and Υ β . The rst one is a linearly spaced grid in the ambient space and Υ β is a regular sampling of the half unit sphere in R 3 . The number of points of the two grids is respectively Λ j

x and Λ β and Λ j is their product. For every nite-dimensional space W * Λ j , it is possible to dene a block matrix K W,Λ j whose blocks are the

RKHS kernel K W,Λ j (x u , ← → β k ), (•, •) u=1,...,Λ j x k=1,...,Λ β between every possible combination of the couples {x u , ← → β k }. The matrix K W,Λ j has a dimension of [Λ j xΛ j ].
We can then use the closest neighbour projection as explained in Sec.2.2.2 to project a varifold on this grid. The squared norm of the projection of the varifold dened in Eq.4.10 results:

|| l δ (p l , ← → n l ) |n l | 2 || 2 W * Λj = u∈Υx p∈Υx k∈Υ β q∈Υ β c uk K W,Λ j (x u , ← → β k ), (x p , ← → β q ) c pq (4.14)
where the scalars c uk and c pq refer to the values resulting from the projection in (x u , ← → β k ) and (x p , ← → β q ) respectively.

Weighted currents The nite dimensional space Q * Λ j for weighted currents is also dened as the span of a set of delta Dirac weighted currents: Span{δ α (y l ,f a m ,f b n ) } where α ∈ R 3 and both the points {y l } and the endpoints {f a m }, {f b n } are constrained to belong to predened grids: y l ∈ Υ j y , f a m ∈ Υ j a and f b n ∈ Υ j b . The three grids are composed of points linearly spaced in the ambient space and their sizes are equal to Λ j y , Λ j a , Λ j b respectively. The product of their sizes is equal to Λ j . As for varifolds, we can dene a block matrix K Q,Λ j for every nite-dimensional space Q * Λ j , whose blocks are the RKHS kernel

K Q,Λ j (y l , f a m , f b n ), (•, •, •)
|| G g=1 w(x g , f a , f b ) T α g || 2 Q * Λj = l∈Υ j y m∈Υ j a n∈Υ j b p∈Υ j y q∈Υ j a r∈Υ j b α lmn K Q,Λ j (y l , f a m , f b n ), (y p , f a q , f b r ) α pqr (4.15)
Chapter 4. A deformation framework to unify morphometry and structural connectivity analysis where the 3D vectors α lmn and α pqr refer to the values resulting from the projection in (y l , f a m , f b n ) and (y p , f a q , f b r ) respectively. The projected norm in Eq.4.15 can be computed using, for instance, the numerical scheme p3m proposed in Durrleman (2010) and based on FFT (Fast Fourier Transform).

Likelihood Once all shapes {S ij } i=1...N and the template T j are projected onto

a nite dimensional space (W * Λ j or Q * Λ j
) it is possible to compute the norm of the projected residuals and the likelihood of the noise. The latter is dened for white and gray matter structures respectively as:

p(ε W ij |σ 2 j ) ∝ 1 |σ 2 j | Λ j /2 exp - 1 2σ 2 j ||S ij -φ All i φ W i (T W j ) || 2 Q * Λj (4.16) p(ε G ij |σ 2 j ) ∝ 1 |σ 2 j | Λ j /2 exp - 1 2σ 2 j ||S ij -φ All i (T G j )|| 2 W *

Λj

where we have scaled the kernels K W,Λ j and K Q,Λ j with the scalar σ 2 j . As it is possible to notice, these equations are very similar to the ones in Eq.4.3 where we modelled the structures as landmarks. The only dierences concern the metric and the denition of the parameter Λ j . These are also the only changes to take into consideration in Eq.4.9 when modelling the structures as varifolds and weighted currents.

In practice, due to the important computational load and time required to cal-

culate both || • || 2 Q * Λj and || • || 2 W * Λj , we choose to compute directly ||(•)|| 2 Q * and ||(•)|| 2 W *
with a fast GPU (CUDA) implementation. Moreover, choosing the number of points for each grid is a tricky task. This becomes even harder when mixing dierent computational models such as varifolds and weighted currents. Here, we propose to select Λ j as the number of points required to build a regular linearly-spaced grid covering the ambient space where the shapes and the template of structure j lie. The distance between every couple of points is equal to λ W or λ g for varifolds and weighted currents respectively. In Sec.4.3.2 we will show that results are stable varying the size of the grids and that we obtain reasonable results with the proposed heuristic rule.

Dieomorphic deformations

We dene here how to compute the dieomorphic deformations of the template complex. Our approach relies on the Large Deformation Dieomorphic Metric Mapping . Calling φ i (x, t) = φ it (x) = x i (t) the position of a point at time t which was located in x at time t = 0, each ow is built by integrating:

∂φ i (x,t) ∂t = v i (φ i (x, t), t) = v i (x i (t), t)
4.2. Methods over t ∈ [0, 1] where v i (x i (t), t) is a time-varying vector eld representing the instantaneous velocity of a point located in x i (t) at time t. Both vector elds v All i and v W i belong to the same RKHS D with Gaussian kernel K D . They are dened by two dierent sets of 3D control points, c All and c W , shared among all subjects, and by two distinct sets of 3D vectors, called momenta, α All i and α W i linked to the control points and specic to each subject i:

v All i (x i (t), t) = K D (x i (t), c All (t))α All i (t) and v W it (x i (t)) = K D (x i (t), c W (t))α W i (t)
, where x i (0) = x and K D (x i (t), c(t)) represents a block matrix of Gaussian kernels with an equal xed bandwidth for both v All i and v W i . The deformation of every point x in the ambient space depends on its initial position at t = 0 and on the evolution of the system L All i (t) = {c All (t), α All i (t)} if the point belongs to the grey matter and on both systems L W i (t) = {c W (t), α W i (t)} and L All i (t) if the point belongs to the white matter. At t = 0 the deformations φ W i0 and φ All i0 are equal to the identity transformations. For both systems, the path from

φ All i0 (resp. φ W i0 ) to φ All i1 (resp. φ W i1
), the latter being the deformation of interest, is chosen as the geodesic one, which means the one that minimizes the total kinetic energy along the path:

1 0 ||v All it || 2 D (resp. 1 0 ||v W it || 2 D ).
ċi (t) = K(c i (t), c i (t))α i (t) = F c (c i (t), α i (t)) c i (0) = c(0) = c 0 αi (t) = -α i (t) T α i (t)∇ 1 K(c i (t), c i (t)) = F α (c i (t), α i (t)) α i (0) = α i0 (4.17)
which can be summarized as

LAll i (t) = F (L All i (t)) (resp. LW i (t) = F (L W i (t))
). The last dieomorphisms φ All i1 and φ W i1 are completely parametrized by the initial conditions of the systems:

L All i (0) = L All i0 = {c All 0 , α All i0 } (resp. L W i (0) = L W i0 = {c W 0 , α W i0 }).
Thus, in order to put into correspondence the template complex T with the subject shape complex S i , we rst integrate forward in time LW i (t) starting from L W i0 and we use these values to deform only the white matter structures of the template complex T W integrating forward in time:

Ṫ W i (t) = K(T W i (t), c W i (t))α W i (t) = G[T W i (t), L W i (t)] T W i (0) = T W i0 = T W (4.18)
The deformed white matter template φ W i1 (T W ) = T W i1 , together with the un- deformed gray matter template T G = T G i0 are then transformed by the second dieomorphism φ All i1 integrating rst forward in time LAll i (t) and then:

Ṫ All i (t) = K(T All i (t), c All i (t))α All i (t) = G[T All i (t), L All i (t)] T All i (0) = T All i0 = T W i1 ∪T G i0 (4.

Optimization procedure

The double dieomorphic framework presented in the previous sub-section can be easily integrated in the Bayesian setting for atlas construction explained in Sec.4.2.2. The two sets of initial control points and momenta {c All 0 , {α All i0 } i=1,...,N } and {c W 0 , {α W i0 } i=1,...,N } represent the deformation parameters to warp the template complex T towards the subject shape complex S i . The initial control points c All 0 and c W 0 are shared among the whole population and we assume they follow a non-informative prior distribution. The initial momenta α All i and α W i are spe- cic to each subject i and they follow a Gaussian distribution exactly as in Eq.4.2.

Modelling the grey matter structures as varifolds and the white matter bundles as weighted currents, assuming that all random variables are independent and considering {α W i0 } i=1,...,N , {α All i0 } i=1,...,N , c All 0 and c W 0 as hidden variables, we can rewrite the cost function in Eq.4.9 as:

M W j=1 N i=1 1 2σ 2 j ||Π(S ij -φ All i1 φ W i1 (T W j ) )|| 2 Q * Λj + P j w j N + M W j=1 1 2 (w j + Λ j N ) log(σ 2 j )+ M G j=1 N i=1 1 2σ 2 j ||Π(S ij -φ All i1 (T G j ))|| 2 W * Λj + P j w j N + M G j=1 1 2 (w j + Λ j N ) log(σ 2 j )+ 1 2 N i=1 (α All i0 ) T (Γ All α ) -1 α All i0 + (w All α + N ) 2 log(|Γ All α |) + w All α 2 tr((Γ All α ) -1 P All α )+ 1 2 N i=1 (α W i0 ) T (Γ W α ) -1 α W i0 + (w W α + N ) 2 log(|Γ W α |) + w W α 2 tr((Γ W α ) -1 P W α ) (4.20)
where M W and M G are the numbers of white matter and grey matter structures respectively. The EM algorithm employed here is an approximation of the Maximum

ΓAll α = N i=1 (α All i0 )(α All i0 ) T + w All α (P All α ) T (w All α + N ) (σ 2 j ) G = N i=1 ||Π(S ij -φ All i1 (T G j ))|| 2 W * Λj + w j P j (w j + N Λ j ) ΓW α = N i=1 (α W i0 )(α W i0 ) T + w W α (P W α ) T (w W α + N ) (σ 2 j ) W = N i=1 ||Π(S ij -φ All i1 φ W i1 (T W j ) )|| 2 Q * Λj + w j P j (w j + N Λ j )
where (σ 2 j ) G and (σ 2 j ) W refer to the grey and white matter structures respectively. The two covariance matrices ΓAll α and ΓW α are equal to a weighted sum between the sample covariance matrix of the initial momenta and the prior. We choose

P All α =K -1 D (c All 0 , c All 0 ) and P W α =K -1 D (c W 0 , c W 0 )
which are block matrices of Gaussian kernels between the initial control points. The matrix K D is the kernel of the RKHS to which belong both vector elds v All i and v W i . This choice is motivated by the fact that if the number of subjects N is denitely smaller than w All α and w W α

then ΓAll α ∝ K -1 D (c All 0 , c All 0 ) and ΓW α ∝ K -1 D (c W 0 , c W 0 )
. This means that the framed regularity terms in Eq.4.20 become

N i=1 (α All i0 ) T K D (c All 0 , c All 0 )α All i0 = N i=1 ||v All i0 || 2 D and N i=1 (α W i0 ) T K D (c W 0 , c W 0 )α W i0 = N i=1 ||v W i0 || 2
D which are the sums of the energies of the geodesic paths of all subjects for both dieomorphisms. This kind of regularity term is often employed in deterministic atlas constructions, namely not based on a statistical setting [START_REF] Joshi | Landmark matching via large deformation dieomorphisms[END_REF][START_REF] Beg | Computing Large Deformation Metric Mappings via Geodesic Flows of Dieomorphisms[END_REF][START_REF] Durrleman | Topology Preserving Atlas Construction from Shape Data without Correspondence Using Sparse Parameters[END_REF]. Furthermore, the use of these priors regularises the estimate of the covariance matrices since it makes possible their inversion even when the number of subjects N is smaller than three times the number of controls points (i.e. size of α All i0 and α W i0 ). This also implies that they can be directly employed in linear statistical analyses such as PCA at the end of the atlas construction without being regularised a posteriori.

The two other parameters, (σ 2 j ) G and (σ 2 j ) W , are equal to a weighted sum between the data-term of the j-th structure and the prior P j . Every parameter balances the importance of structure j with respect to the other structures and with respect to the regularity terms of both dieomorphisms. The chosen prior imposes a minimum value to σ2 j which is useful to avoid overtting. In fact, without a prior, the min- imisation process might focus only on a structure k reducing its residuals almost to zero and ignoring the other structures. This would result in σ2 k → 0 and therefore also in log(σ 2 k ) →-∞. The gradients of the cost function E in Eq.4.20 with respect to the other parameters are equal to: Chapter 4. A deformation framework to unify morphometry and structural connectivity analysis

∇ T G k E = N i=1 1 2σ 2 k ∇ T G k D G ik ∇ T W k E = N i=1 1 2σ 2 k ∇ T W k D W ik ∇ α All s0 E = M j=1 1 2σ 2 j ∇ α All s0 D sj + (Γ All α ) -1 α All s0 ∇ α W s0 E = M j=1 1 2σ 2 j ∇ α W s0 D sj + (Γ W α ) -1 α W s0 ∇ c All 0 E = N i=1 M j=1 1 2σ 2 j ∇ c All 0 D ij ∇ c W 0 E = N i=1 M j=1 1 2σ 2 j ∇ c W 0 D ij (4.21)
where

D G ij = ||Π(S ij -φ All i1 (T G j ))|| 2 W * Λj and D W ij = ||Π(S ij -φ All i1 φ W i1 (T W j ) )|| 2 Q * Λj
refer to the data terms of the grey and white matter structures respectively whereas D ij refers to the data-term of any structure. Moreover, we have not taken into account the dierentiation of the prior P α with respect to the control points since its norm is negligible with respect to the one of the data term.

To calculate the gradients of the data terms {D ij } we need to compute the deformed template complex obtained at the end of the second dieomorphism for every subject i and for all structures, namely φ All i1 (T G ) and φ All i1 φ W i1 (T W ) . First, we integrate forward in time LW i (t) = F (L W i (t)) (Eq.4.17) and

Ṫ W i (t) = G[T W i (t), L W i (t)]
(Eq.4.18). Then, we integrate LAll i (t) = F (L All i (t)) (Eq.4.17) and, using as initial value T All i (0) = T W i1 ∪ T G i0 , we integrate also 4.19). After that, we can compute the data term D ij and its gradient with respect to the vertices of T All i (1) = T All i1 . The computations of the gradient for surfaces mod- elled as varifolds are shown in the Appendix of [START_REF] Durrleman | Morphometry of anatomical shape complexes with dense deformations and sparse parameters[END_REF] whereas the calculations about the gradient for ber bundles modelled as weighted prototypes are presented in Sec.4.A.1. Using the calculus of variations (see Sec.4.A.2), the information about the gradient is brought back from t = 1 to t = 0 to update rst L All i (0) = {c All 0 , α All i0 } and T G and then L W i (0) = {c W 0 , α W i0 } and T W . The optimisation is based on a set of linearised ODEs describing the evolution of four Instead, the size of θ W i and ξ W i are equal to the ones of T W and L W i respectively.

Ṫ All i (t) = G[T All i (t), L All i (t)] (Eq.
auxiliary variables θ All i , ξ All i = {ξ All αi , ξ All ci }, θ W i , ξ W i = {ξ W αi , ξ W ci }: θAll i (t) = -(∂ T All i G All i (t)) T θ All i (t) θ All i (1) = ∇ T All i (1) D i (4.22) ξAll i (t) = -(∂ L All i G All i (t)) T θ All i (t) + (d L All i F All i (t)) T ξ All i (t) ξ All i (1) = 0 (4.23) θW i (t) = -(∂ T W i G W i (t)) T θ W i (t) θ W i (1) = θ All,W i (0) (4.24) ξW i (t) = -(∂ L W i G W i (t)) T θ W i (t) + (d L W i F W i (t)) T ξ W i (t) ξ W i (1) = 0
We rst integrate backward Eq.4.22 and Eq.4.23 obtaining θ All i (0) and ξ All i (0) = {ξ All αi (0), ξ All ci (0)}. Then, we use θ All,W i (0), which are the initial values of θ All i relative to the white matter structures, as nal values for θ W i and we integrate backward Eq.4.24 and Eq.4.25 obtaining θ W i (0) and ξ W i (0) = {ξ W αi (0), ξ W ci (0)}. From this set of equations, we can notice that the optimisation of the two dieomorphisms is linked by the constraint θ W i (1) = θ All,W i (0) and that the information given by ∇ T All i (1) D i where D i = {D ij } j=1,...,M ows rst through the second dieomorphism All and then through the rst one W and eventually it is used to update all the parameters in the following way: For surfaces, we use the average of the population when it is available a vertexcorrespondence, otherwise we use a regular sphere which is rst centred with respect to the ensemble of subject meshes and then it is scaled to an ellipsoid using the three main modes of variations of the meshes considered as points cloud.

∇ T G E = N i=1 θ All,G i (0) ∇ T W E = N i=1 θ W i (0) ∇ α All s0 E = ξ All αs (0) + (Γ All α ) -1 α All s0 ∇ α W s0 E = ξ W αs (0) + (Γ W α ) -1 α W s0 ∇ c All 0 E = N i=1 ξ All ci (0) ∇ c W 0 E = N i=1 ξ W ci (0)
For bundles of weighted prototypes, we propose to gather all the original streamlines of the subjects into a single bundle which is then approximated as a set of weighted prototypes. The weights of the prototypes are then scaled so that the norm of the template is equal to the average norm of the population.

Experiments and Results

The dataset used in this chapter contains 76 subjects: 27 controls and 49 patients subject to Gilles de la Tourette syndrome (GTS) and divided in three sub-groups based on their symptoms: G1=simple-tics, G2=complex-tics, G3=complex tics with obsessive compulsive disorders (OCD).

In this section, we will rst evaluate the robustness of the proposed framework with respect to the size of both the varifold and weighted current grid using left caudate and the bundle of prototypes connecting it to the left hemisphere of the cortical surface. Then, we will take advantage of the proposed atlas construction procedure to investigate morphological and organisational abnormalities of the left cortico-putamen circuit in the population of GTS patients. These alterations are likely to be associated to a pathological development of the brain. The employed neural circuit is composed of left putamen, left hemisphere of the cortical surface and the streamlines (prototypes) connecting them.

In the following experiments, we rst compute an atlas of 50 iterations using only the sub-cortical structures employed. This permits to transform the ellipsoid in a structure more similar to the sub-cortical nucleus. After that, we cut the prototypes of the template bundle at the intersection with the inter-linked grey matter structures.

Experiments and Results

Materials

Anatomical scans of all subjects were acquired using 3D T1-weighted sequences with a voxel size of 1x1x1 mm 3 . Diusion weighted scans were acquired with sequences of 50 directions with a B-factor of 1000 and a voxel size of 2x2x2 mm 3 . For each subject, we consider the left hemisphere of the cortical surface, left putamen, left caudate and the ber bundles connecting these sub-cortical nuclei to the left hemisphere of the cortex. The cortical surface is segmented using FreeSurfer v5.3 [START_REF] Fischl | Automatically Parcellating the Human[END_REF]) from the T1-weighted images followed by a pipeline of BrainVisa 4.

Robustness with respect to the grid size Λ j

We evaluate here the robustness of the proposed method with respect to the number of points of the grids (Λ j ) where both the structures and the template modelled as varifolds or weighted currents are projected to. We estimate 16 atlases using a group of 8 anatomical complexes composed of left caudate and left caudate bundle. The left caudate is modelled as varifold using λ W =3mm and the left caudate bundle as weighted currents with λ g =7mm, λ a =10mm (cortex) and λ b =5mm (caudate).

Atlases share the same parameters except for the grid size of both structures. In It is possible to notice that when the grid size of the bundle is greater than Chapter 4. A deformation framework to unify morphometry and structural connectivity analysis the one of the nucleus (W C > V ), results are very stable and satisfactory. The only exception is when the dierence is too great, namely with V = 1000 and W C = 500.000. In such a case, the residuals of the caudate are higher than in the other examples and the template does not evolve since its data-term is too small with respect to the one of the bundle (see Fig. 4.5 top-left). On the contrary, when

V > W C like with V = 100.000 and W C = 1.000, we have the opposite eect since the residuals of the bundle are higher and its template does not converge (see Fig. 4.5 top-right).

In the following experiments, we propose to use an heuristic to dene V and W C. They are selected as the number of points required to build a regular linearlyspaced grid covering only the subject shapes and the template of each structure and the distance between every couple of points is given by the smallest value between λ W and λ g among all structures. In this example, the smallest value is the one of the caudate, namely 3mm, and the grid sizes result V ≈ 5.000 and W C ≈ 110.000.

Results based on this choice are satisfactory and as it is possible to see from Fig. 4.5 at the bottom, both templates are updated during the atlas construction.

Experiments and Results

Group dierences

In all the following experiments we use the cortico-putamen neural circuit. The cortical surface is modelled with landmarks, the putamen as varifold with λ W =3mm and the ber bundle connecting them as weighted currents with λ g =7mm, λ a =10mm (cortex) and λ b =5mm (putamen). All bundles are approximated as weighted prototypes. The bandwidths of both dieomorphic kernels is equal to 11mm and the number of control points is 816 for both dieomorphisms. The maximum number of iterations is 120 and all the computations are performed on a Intel Xeon, 32 cores, CPU E5-2650, 2.60GHz with a graphic card NVIDIA Quadro 5000. The code is written in C++ and CUDA and it is an extension of the freely available software deformetrica (www.deformetrica.org). The computational time for an atlas of 30 subjects is equal to 144 hours for 120 iterations. All shape complexes are previously rigidly registered to a reference shape complex.

It is important to notice that a bandwidth of 11mm permits to deform all white matter ber bundles in a correlated way with at least one grey matter mesh. This permits to have an identiable model with a unique decomposition of the two diffeomorphisms all over the ambient space.

Atlas with pooled subjects

Toy example We start this section with an explanatory example based on a toy data-set constituted of 6 pseudo shape complexes belonging to two dierent populations (i.e. controls and GTS patients). Each complex is composed of a pseudo cortical surface, a pseudo sub-cortical structure and a pseudo ber bundle linking them. They are shown in Fig. 4.6 where it is possible to notice that the complexes of population A have a dier organisation and shape with respect to the ones of population B. Indeed, subjects from A have a bigger cortex and their sub-cortical structure have a bump on the right whereas in population B it is on the left and the ber bundles of the two populations integrate completely dierent grey matter areas. In population A bers go from the right area of the sub-cortical structure to the left area of the cortex whereas in population B it is exactly the contrary.

We estimate an atlas considering both populations together as if they were from a single group. We start with a simple and unbiased initial template complex and we obtain a nal template complex and the covariance matrices of the momenta of both dieomorphisms. The template shows the features common to both populations. The two covariance matrices describe the morphological and organisational variability within the 6 subjects. We compute a Principal Component Analysis (PCA) for each covariance matrix and we deform the nal template complex at ±σ (standard deviation) along the rst modes of both PCAs.

Since the intra-group variations are denitely smaller than the inter-group ones, the main variations among the 6 subjects captured by the white dieomorphisms {φ W i }, which aect only the ber bundles, explain the principal dierences in structural connectivity between the two populations. The positions of the ber bundle Chapter 4. A deformation framework to unify morphometry and structural connectivity analysis at -σ and +σ are the ones of population A and B respectively. Moreover, even the rst mode of the global PCA describes the main global morphological variations between the two populations. As previously, we can notice that both grey matter structures at -σ and +σ reproduce the morphological characteristics of population A and B respectively.

This example shows the exploratory potential of the proposed algorithm even if it is based on a simple toy data-set where the inter-group dierences are probably exaggerate compared to a real-data example. Nevertheless, given the important structural changes that are likely to occur in syndromes such as GTS, we may assume that controls and patients create separate clusters. The rst modes of variations of the two dieomorphisms should therefore describe the main dierences between these two clusters highlighting the eects of the pathology on both the anatomy and structural connectivity.

This could not be achieved using a single dieomorphic atlas since the templateto-subject registrations would not be able to put into correspondence all the structures of the complex. An example of registration is shown in Fig. 4.7, where it is possible to notice that the matching using a single dieomorphism is not accurate.

Instead, a double dieomorphic registration makes rst the ber bundle sliding onto the grey matter structures and then it accurately aligns with the second (global) dieomorphism all structures of the template complex to the homologous ones of the subject complex. This permits to disentangle the dierences in structural connectivity, captured by the rst dieomorphism, from the global morphological changes, captured by the second dieomorphism.

Cortico-putamen circuit Following the same strategy we also build an atlas with pooled subjects based on real data. We use the cortico-putamen circuits of 10 controls and 10 patients subject to GTS. In Fig. 4.8 we show the estimated template complex and its deformations along the rst mode at ±3 std of a PCA based on the initial momenta of the white dieomorphisms. The rst and last column present the densities of the endpoints of the template bundle onto the xed grey matter structures. Colours in the middle column highlight the displacement of the deformed template bundle along the rst mode. The two arrows point to the most variable areas. The rst one, at the top, is in the superior-frontal area of the cortical surface.

As it is possible to notice from the densities on the cortex, bers tend to concentrate in a specic area in the rst row which is then left almost empty in the bottom row.

The second arrow points to the lateral part of the putamen where bers move from the bottom to the central part in the rst and last row respectively. This movement is clearly visible in the right column.

In Fig. 4.9 we present the main global morphological variations of both the cortical surface and the putamen. We compute the rst two modes of a PCA based on the initial momenta of the global dieomorphisms. Then, we deform the estimated template (shown in the middle in gray) along the two modes at ±3 std. Colours highlight the displacement of the deformed template cortex. The rst mode captures changes mainly in the parietal area of the cortex whereas the second mode The changes presented in these two gures could highlight the main variations between the two groups or they could show their common main variability. In the next sub-section we will compare these results with the ones obtained by analysing Chapter 4. A deformation framework to unify morphometry and structural connectivity analysis each group independently.

We also computed the average of the initial momenta of each population and we deformed the template complex along these two directions. This gave us a quantitative way to compare the average congurations of both groups. Results

(not shown) were quite similar and no dierences stood out.

Experiments and Results

Separated Atlases

We compute here two dierent atlases using 30 patients (10G1 + 10G2 +10G3) and 27 controls respectively. We use the same initial template complex as in the previous sub-section. We compute the rst three modes of a PCA based on the initial momenta of the white dieomorphisms for both groups. Then, we deform the nal template complex along every mode at ±3 std. Deformations are shown in Fig. 4.10 for controls and in Fig. 4.11 for patients. These two gures present the main variability in structural connectivity within each group.

By comparing them, we can notice that the rst mode of each group captures variations relative to the commisural bers passing through the corpus callosum.

This means that variations in this area are probably common to both groups. The other two modes show instead changes aecting dierent areas. In particular, in controls the most variable area is the superior-frontal instead in patients the variations are more heterogeneous. Interestingly, the changes of the second mode in patients regard almost the same areas as in Fig. 4.8. This might signify that the group of patients, due to the vast clinical heterogeneity of GTS, might create a large cluster in the shape space which might even intersect the one of controls creating a kind of continuum from the group of controls to the most severe patients. If we suppose then that the intra-variability of the group of patients is greater than the one of controls we could explain why the main variations found considering controls and patients together are very similar to the greatest changes estimated among only patients.

Chapter 4. A deformation framework to unify morphometry and structural connectivity analysis First of all, we build an atlas with 10 subjects (5 controls and 5 patients) which results in a template, two sets of control points and two well-conditioned covariance matrices of deformation parameters of the white (Γ W α ) and global (Γ All α ) dieomor- phism respectively. Since we use subjects from both clinical groups, the template should be positioned in between them in the shape space. The template is successively warped to all the remaining J subjects xing the control points equal to the estimated ones and using the two covariance matrices as regularity terms, namely minimising for each subject:

M W j=1 1 2σ 2 j ||Π(S j -φ All 1 φ W 1 (T W j ) )|| 2 Q * Λj + (α All 0 ) T (Γ All α ) -1 α All 0 M G j=1 1 2σ 2 j ||Π(S j -φ All 1 (T G j ))|| 2 W * Λj + (α W 0 ) T (Γ W α ) -1 α W 0 (4.27)
which is equal to Eq.4.20 but with only one subject (index i disappears) and with a constant σ 2 j which is xed to the value estimated in the atlas. The resulting initial momenta {α All i0 } 1...J and {α W i0 } 1...J represent the features used as input in the classier. We employ Linear Discriminant Analysis (LDA) with a leave-one-out cross validation strategy. We test the discriminative power of the two dieomorphisms separately by using either only {α All i0 } 1...J or only {α W i0 } 1...J . In each experiment, we suppose that the class-conditional densities of the initial momenta for both clinical groups are Gaussian, namely f c (α i0 ) = N (µ c , Γ c ) and f p (α i0 ) = N (µ p , Γ p ), and that they have priors equal to their relative number, namely π c = J c /J and π p = J p /J where the indices c and p refer to controls and patients respectively, π c and π p are the prior distributions, J c = 22 is the number of controls and J p = 44 is the number of patients.

In LDA the feature space is divided into two parts by a decision boundary which is based on the Bayes factor B dened as:

B(α t ) = log f c (α t ) f p (α t ) + log π c π p = - 1 2 (α t -µ c ) T Γ -1 c (α t -µ c ) - 1 2 log |Γ c | + 1 2 (α t -µ p ) T Γ -1 p (α t -µ p ) + 1 2 log |Γ p | + log( J c J ) -log( J p J ) (4.28)
where α t is the initial momenta of the test subject t and we omitted the index

Experiments and Results

0 for clarity purpose. If B is smaller than a threshold τ the subject is classied as patient otherwise as control. Moreover, in LDA the covariance matrices of the two groups are assumed to be equal and with full rank (i.e. Γ c = Γ p = Γ) which means that Eq.4.28 becomes:

B(α t ) = α T t Γ -1 (µ c -µ p ) w - 1 2 µ T c Γ -1 µ c - 1 2 µ T p Γ -1 µ p -log( J c J ) + log( J p J ) b (4.29)
The decision boundary is therefore linear, namely it is the hyperplane in the feature space: α T t w -b -τ = 0. Depending on the side where the α t lies, the correspondent test subject t is classied as patient < 0 or control > 0.

To estimate the accuracy of our model without incurring in over-tting, we use a leave-one-out cross validations technique. At each iteration, we estimate the best hyperplane using J-1 subjects and we test it on the left out subject t. First, we compute the averages of both controls and patients within the J-1 subjects. Assuming that the test subject t is a control, it results:

µ c = 1 J c -1 Jc-1 i=1 α ic µ p = 1 J p Jp i=1 α ip (4.30)
After that, using the covariance matrix estimated in the initial atlas as Γ in Eq.4.29, we compute the best threshold τ as the one which gives the highest balanced accuracy BA, namely the average between specicity and sensitivity. When there is a continuous range of values at which BA is maximal, we choose the middle one.

Once estimated the best hyperplane in the training set, we use it to test the left-out subject t. We iterate this process for all the J subjects of the data-set. Resulting sensibility, sensitivity and balanced accuracy are shown in Table 4.1.

Since the group of patients is composed of three sub-groups, namely G1=simpletics (17 patients), G2=complex tics (15) and G3=OCD (12), we test our classication algorithm using these sub-groups separately or together. Moreover, we compare these results with the ones obtained using the initial momenta of a single dieomrphism. Following the same strategy, we rst estimate a single-dieomorphic atlas (like in Chapter 2) using the same 10 subjects as before. Then, the estimated template is warped with a single-diemorphic deformation towards the remaining J subjects. Even in this case we x the control points and we use the estimated covariance matrix as regularity term. Eventually, we employ LDA with a leave-one-out cross validation technique to classify the initial momenta exactly as before.

As it is possible to notice from Table 4.1, the classication scores based on a single dieomorphism are worse than using the proposed double dieomorphism especially when using the most severe groups of patients (G2 and G3).

Due to the variability of the results in Table 4.1, we also investigate their distributions within the group of patients with a bootstrap analysis. More precisely, we Chapter 4. A deformation framework to unify morphometry and structural connectivity analysis perform it on the top of the previous leave-one-out cross validation classication.

At each iteration of the bootstrap, we pick a random sample (with replacement) of the 44 patients which is classied, together with the (un-modied) 22 controls, using LDA. We repeat this process 1000 times. Each bootstrap sample gives as result an estimate of sensitivity, sensibility and balanced accuracy. The histograms of balanced accuracy for both the double dieomorphism and the single dieomorphism are shown in Fig. 4.12. The averages for sensitivity and specicity are respectively: 74% and 51% for the global dieomorphism and 73% and 64% for the white dieomorphism. Instead, the averages sensitivity and specicity for the single dieomorphism are 64% and 48% respectively.

The results obtained with the double dieomorphism are, even in this case, better than using a single dieomorphism. This demonstrates that the proposed double dieomorphic approach captures important and relevant information which better characterise the pathological anomalies of GTS. Moreover, looking at Table 4.1, we can see that there is a general trend in the classication scores. Whatever the dieomorphic scheme employed, we always classify the most severe patients, namely groups G2 and G3, better than the ones with only simple tics (G1). This conrms the hypothesis that there might be a continuum in the shape space between healthy controls and GTS patients where the subgroup G1 would be closer to the healthy group than the subgroup G3. Furthermore, it also conrms that the atlas construction technique might be a good model to investigate GTS and more generally neurodevelopmental syndromes.

Discussion and Conclusion

We presented a method to analyse the organisational and morphological variability of a population of anatomical complexes composed of grey matter surfaces interconnected by white matter ber bundles. It is built on a double-dieomorphic mesh-based atlas construction embedded in a Bayesian framework. The resulting atlas consists of a template complex, showing the common characteristics of the population, and two well-conditioned covariance matrices describing the distributions of the deformation parameters, namely the initial momenta, of both dieomorphisms.

The cascade of two dieomorphisms permits to put into correspondence anatomical complexes characterised by ber bundles connecting dierent areas of the grey matter, namely with a dierent structural connectivity. The rst (white) dieomorphism modies the relative position of the ber bundles with respect to the xed grey matter structures in such a way that it is possible to correctly align the entire brain, both grey and white matter structures, with the second (global) dieomorphism. In this way, we can disentangle the morphological variations into two categories: the ones related to the structural connectivity are captured by the rst dieomorphism and the ones about global morphological dierences, aecting both white and grey matter structures, are instead described by the second dieomorphism. The deformations obtained with the white dieomorphisms can be compared across subjects since they are computed with respect to the same reference frame, that it the xed grey matter. On the contrary, it should be noted that the deformations of the ber bundles obtained with the global dieomorphisms can not be compared across subjects. Global dieomorphisms can be used to study the morphological variability of the grey matter structures. Both white and global changes can be visualised as a continuous deformation of the template complex. This permits to easily characterise, localise and quantify both organisational and morphological pathological anomalies altering grey and white matter structures.

Furthermore, our algorithm is limited for now to cases in which we can assume that there is a point-correspondence among the cortical surfaces in the population.

In future works, we will try to adapt our technique for considering also subjects showing dierent cortical gyrications. Another possible extension would be the use of a probabilistic tractography instead than the deterministic one employed here. This would make the analysis computationally more demanding but thanks to the approximation scheme presented in Chapter 3 it should still be feasible. ). These methods can be considered as data-driven and they seem to work well in practice. However, due to the importance of the initialisation in gradientdescent algorithms, we think that it would be relevant to investigate the robustness of the results with respect to a change in the initial template. For instance, we could randomly divide a group of subjects into dierent sub-groups. Then, we could create the initial template and estimate the atlas in each sub-group separately. Results could then be compared across sub-groups giving an estimate of the robustness.

However, this would require a considerable computational time.

All experiments shown in this chapter were about one neural circuit. To full the goal of this thesis, namely an holistic approach integrating the whole CSTC network in a single framework, we should use the proposed method with all the neural circuits at the same time. However, this might not work well since every circuit might be aected by dierent pathological alterations. This means that every bundle should be deformed in a dierent and independent way with respect to the others. Thus, we should not use a single dieomorphism to deform the entire white matter when composed by multiple bundles. A possible solution would be to use N +1 dieomorphisms. Every bundle would be independently deformed by a white dieomorphism. Then, all bundles would be deformed, together with the grey matter structures, by another global dieomorphism. In this way, we could capture the dierent variations in structural connectivity proper to each bundle and also the global morphological changes associated to all structures of the CSTC circuits.

In this chapter, we used three dierent computational models: landmarks, varifolds and weighted currents. When using the last two innite-dimensional spaces,

Discussion and Conclusion

we need to project both the subject shapes and the template to a nite-dimensional space where probability density functions can be computed. We assumed that the noise of each structure j follows a Gaussian distribution whose covariance matrix has the size of the nite-dimensional approximation space for weighted currents and varifolds and it is equal to the number of point-correspondences when working with landmarks. In both cases we dene it as Λ j .

Using ω j =0.01Λ j N and P j = 0.05R 0 j ω j with R 0 j equal to the initial residual of struc- ture j, as proposed in Sec.2.3.2, we can notice that the data-terms of the cost function of the atlas construction (Eq.4.20) can be approximated by:

≈ M j=1 N Λ j 2 log R j + 0.05R 0 j N Λ j (4.31)
after removing the constant terms, with R j equal to the residual of structure j and where N and M are the number of subjects and structures respectively. This shows that the data-term of each structure j depends on the metric of the computational model chosen and on the value of Λ j . If all structures are modelled with the same computational model, like in Chapter 2, the data-terms are commensurable since we use the same rule to compute both R j and Λ j but if the structures are modelled with dierent computational models, like in this chapter, nothing guarantees that the data-terms are commensurable. Depending on the number of landmarks chosen or on the size of the grid where a varifold or weighted current is projected to, a structure could be over-weighted with respect to the others (or the contrary). In this chapter, we used some heuristics to nd a good balance between the dierent data-terms which gave reasonable results but future works will have to nd a way to harmonize the choice of Λ j among the dierent computational models.

In the proposed method, we assumed that the initial momenta of the two dieomorphisms are independent, that is to say that p(α All i , α W i ) = p(α All i )p(α W i ), even if the update rule for α All i and α W i are related as explained in Sec.4.2.5. It would seem more reasonable to take that into account by modelling directly p(α All i , α W i ) without the assumption of independence. We could model, for instance, their joint distribution as a single Gaussian distribution. However, the statistical relationship between α All i and α W i is highly complex since they are related by the linearised ODEs shown in Sec.4.2.5 and we have not found yet a satisfactory solution to model their joint distribution. This is left as future work.

Nevertheless, we demonstrated that the proposed double dieomorphic approach can detect some useful and relevant information since it leads to better classication scores than a single dieomorphism. This means that the information about structural connectivity might play an important role in the characterisation of GTS and that the proposed method can bring an important contribution in the description of the pathophysiological mechanisms underlying GTS. Q * and then its gradient with respect to every point of the ber F k = L k n=1 δ α kn (x kn ,f k ,t k ) where x kn is the center of segment n of ber k and α kn is its tangent vector. They are dened as x kn = (p kn + m kn )/2 and α kn = p kn -m kn , where p kn and m kn are the rst and the last point of segment n of ber k respectively. This means that every couple of neighbour segments n and n + 1 has one point in common: p n = m n+1 (except for the extremities). We assume also that all bers are oriented from the extremity f i (v p ) towards t i (u p ) and we call k a the kernel related to the rst extremity (f i ) and k b the kernel for the second extremity (t i ). Since weighted currents are dened with respect to the center and tangent vectors, we need to use the chain rule. Moreover, every point (except for the extremities) depends on two segments since it can be seen as p n but also as m n+1 which means that its update rule will depend on the center and tangent vectors of both segments n and n + 1. The derivative with respect to the point p ks (resp. m ks ), that is

∂|| N i=1 τ i F i -M p=1 wpTp|| 2 Q * ∂p ks (resp. ∂||•|| 2 Q * ∂m ks
), if it is not the extremity (p ks = t k , resp.

2τ k N i=1 i =k τ i k a (f i , f k )( -2(t k -t i ) λ 2 b )k b (t i , t k ) F k , F i U * - 2τ k M p=1 w p k a (v p , f k )( -2(t k -u p ) λ 2 b )k b (t k , u p ) F k , T p U *
Instead when computing the derivative with respect to m ks and m ks = f k one needs to add:

2τ k N i=1 i =k τ i k a (f i , f k )( -2(f k -f i ) λ 2 a )k b (t i , t k ) F k , F i U * - 2τ k M p=1 w p k a (v p , f k )( -2(f k -v p ) λ 2

4.A.2 Gradient of the Atlas construction procedure

We compute here the gradients of the atlas construction as explained in Sec.4.2.5.

A variation of the initial control points and momenta of the rst dieomorphism δL W i0 produces a variation in their evolution δL W i (t) and consequently a variation of the path of the ow of the rst dieomorphism δT W i (t). In parallel, δL All i0 produces a variation in δL All i (t) which , together with δT W i1 , induces a variation in the path of the ow of the second dieomorphism δT All i (t) and consequently also in the criterion δE.

δ LW i (t) = (d L W it F W i (t)) T δL W it δL W i (0) = δL W i0 δ Ṫ W i (t) = (∂ T W it G W i (t)) T δT W it + (∂ L W it G W i (t)) T δL W it δT W (0) = δT W 0 δ Ṫ All i (t) = (∂ T All it G All i (t)) T δT All it + (∂ L All it G All i (t)]) T δL All it δT All i (0) = δT W i1 ∪ δT G 0
The two rst ODEs are linear and their solution is given by: δL(t) = exp

Conclusive summary

In this thesis, we proposed a new approach to localise, characterise and visualise the anatomical abnormalities associated to a neurodevelopental disorder and aecting both the structural organisation and shape of neural circuits. This methodology is based on the Grenander's atlas construction procedure which permits to analyse the neural circuits as a single interwoven complex where every component is virtually represented as a 3D mesh. This makes possible the integration of white matter ber bundles and grey matter structures in a single framework. Moreover, the use of meshes facilitates the visualisation of the complex and the interpretation of its pathological alterations.

In Chapter 2 we proposed a Bayesian framework in which embed the generative model of the atlas construction procedure. It allows to automatically estimate the noise variance of each component of the neural circuits which should otherwise be xed by the user. We showed that these values can highly inuence the results and their automatic estimate becomes even more crucial if the number of analysed structures increases. This is therefore essential to perform an holistic analysis of the entire neural circuits. Moreover, the proposed statistical setting also permits to estimate a well-conditioned covariance matrix of the deformation parameters. This makes thus possible to directly use it for statistical analysis, such as PCA, without regularising it a posteriori.

Furthermore, we proposed to model every surface of the grey matter as a Gaussian random varifold for which we dene nite-dimensional approximation spaces similarly to [START_REF] Bibliography Durrleman | Statistical models of currents for measuring the variability of anatomical curves, surfaces and their evolution[END_REF] where it is possible to dene probability density functions. This computational model does not need point-correspondences and it has a closed-form metric which is easily derivable. Moreover, it can be used for any kind of shape, including 3D curves such as the white matter bers. However, this framework takes into consideration only the geometry (pathway) of the bers. This makes the template-to-subject registrations accurate only in the most dense parts of the bundles and it almost ignores the extremities, namely where the streamlines integrate grey matter structures. Thus, it hampers the study of variations in structural connectivity. In addition to that, ber bundles resulting from tractography algorithms would require an enormous computational load and time to be directly used in the atlas construction procedure due to their huge number of streamlines.

To this end, we proposed in Chapter 3 an approximation scheme for ber bundles.

The resulting parsimonious and computationally-ecient representation makes feasible the integration of ber bundles in the atlas procedure together with the other components of the neural circuits. Moreover, it is visualised as a scattered set of tubes whose radii are proportional to the number of streamlines approximated. The smaller number of tubes, with respect to the original amount of streamlines, eases the visualisation and interpretation of the bundles since it allows to clearly identify and analyse all its fascicles. Moreover, it drastically reduces the computational load for visualisation software thus accelerating any kind of manipulation. The proposed approximation is based on the computational model of weighted currents, an adaptation of the one of functional currents (Charon and Trouvé, 2013a). Similarly to currents, it does not need point correspondences or streamlines correspondences and it augments its denition taking into consideration not only the pathway of the streamlines but also the locations of their extremities. This makes thus possible to correctly register also the extremities of the bundles in the template-to-subjects deformations and, theoretically, to capture the variations in structural connectivity.

However, the deformations of the template complex in standard atlas construction procedures are dened as single dieomorphisms of the entire ambient space.

This precludes the study of variations in topology such as the changes in structural connectivity. To this end, we proposed in Chapter The proposed method describes the organisation and morphology of a group of anatomical complexes composed of both grey matter surfaces and white matter tracts. We showed its eectiveness by comparing the cortico-putamen circuits of a group of patients subject to Gilles de la Tourette syndrome and a control group.

Limitations and Perspectives

Preliminary results highlight abnormal dierences about both the shape of the grey matter structures and the structural connectivity. These alterations are in line with the hypothesis put forth in [START_REF] Worbe | Distinct structural changes underpin clinical phenotypes in patients with Gilles de la Tourette syndrome[END_REF][START_REF] Worbe | Altered structural connectivity of cortico-striatopallido-thalamic networks in Gilles de la Tourette syndrome[END_REF]. Moreover, we also showed that the proposed approach leads to better classication scores than the single dieomorphic method described in Chapter 2. This suggests that it might better characterise the anatomical alterations associated to GTS and therefore that it could bring an important contribution in the description of the pathophysiological mechanisms underlying GTS.

The algorithms presented in this thesis have been coded mostly in C++ and some parts also in CUDA. The code of the proposed Bayesian framework (Chapter 2) and of the double dieomorphic deformation (Chapter 4) will be integrated to the software deformetrica www.deformetrica.org.

Limitations and Perspectives

The proposed methodology oers a mathematical tool which can help understanding the pathophysiology of neurodevelopmental syndromes. This opens a wide spectrum of possible future clinical applications such as autism and attention decithyperactivity disorders.

Moreover, to fully characterise the vast clinical heterogeneity and complexity of neurodevelopmental disorders we still need to extend the proposed method. This oers a challenging opportunity for further methodological improvements. In the following, we present perspectives and possible solutions to some of the limitations highlighted in the previous chapters.

Mixture of templates

A possible extension would be to use a mixture of templates instead than a single one as proposed in [START_REF] Allassonnière | Towards a Coherent Statistical Framework for Dense Deformable Template Estimation[END_REF]. Statistical learning techniques could be used to automatically estimate the number of templates. Each template would characterise a sub-group of the population. This could then be used to cluster the patients in smaller groups which would be (hopefully) clinically more homogeneous.

The analysis of the templates and their variability could help understanding whether these sub-groups share anatomical alterations, distinctive therefore of the disorder, or if they are characterised by specic abnormalities. Eventually, this could help clarifying the pathophysiological mechanisms related to the dierent symptoms.

Cortical surface

In this thesis, we modelled the cortical surface with landmarks which is a computationally convenient solution but it is based on a strong hypothesis of pointcorrespondence which might not be correct. A dierent strategy would be modelling it as a varifold with explicit sulcal constraints given, for instance, by the 120 Chapter 5. Conclusions and Perspectives BrainVisa sulci [START_REF] Perrot | Cortical sulci recognition and spatial normalization[END_REF][START_REF] Auzias | Dieomorphic Brain Registration Under Exhaustive Sulcal Constraints[END_REF], the 56 major sulci of MNI Display [START_REF] Holmes | Cortical peeling: CSF/grey/white matter boundaries visualized by nesting isosurfaces[END_REF][START_REF] Joshi | Dieomorphic Sulcal Shape Analysis on the Cortex[END_REF] or the so called sulcal pits, namely the deepest points of the folding pattern [START_REF] Auzias | Deep sulcal landmarks: Algorithmic and conceptual improvements in the denition and extraction of sulcal pits[END_REF]. In this way, we would account for the cortical shape variation and at the same time we would be able to align the cortical folding architecture.

However, we should also take into consideration the fact that patients with Gilles de la Tourette syndrome might have an abnormal structural pattern of cortical sulci, as it was demonstrated in [START_REF] Muellner | Altered structure of cortical sulci in gilles de la Tourette syndrome: Further support for abnormal brain development[END_REF].

Integrating functional information

Another possible extension would be to perform an even more holistic approach by adding a functional signal to the anatomical complexes. This could be easily taken into account by augmenting the framework of varifolds and weighted currents as in functional currents (Charon and Trouvé, 2013a). The functional signal could be represented, for instance, by functional MRI or EEG. The functional information would be directly taken into account by augmenting the dimensionality of the space where shapes lie. It would change from a shape space to a functional-shape space where the deformation of the template complex towards a subject complex would be driven by a similarity measure which considers both shape and functional signal.

This approach would be in contrast with most of the works present in the literature which either study functional and structural imaging separately or use the latter to restrain the functional analysis to specic areas [START_REF] Schumann | Abnormal structure or function of the amygdala is a common component of neurodevelopmental disorders[END_REF][START_REF] Vaidya | Neurodevelopmental Abnormalities in ADHD[END_REF][START_REF] Singh | Functional and structural abnormalities associated with Bibliography empathy in patients with schizophrenia: An fMRI and VBM study[END_REF].

The framework of functional currents could also be used to augment the approximation scheme proposed in Chapter 3 with information about the white matter microstructure. Every streamline could be characterised by a vector of scalar quantities, such as FA or GFA, describing the white matter structural integrity of the voxels through which the ber has been traced. A similarity measure between two streamlines would then take into account their pathways, connectivity and also the local white matter microstructure. This could be combined with a new deformation setting always in the spirit of the Grenander's approach like in [START_REF] Charlier | The fshape framework for the variability analysis of functional shapes[END_REF].

The dierences between two streamlines would always be described by warping the rst one onto the other. However, in this case the deformation would change not only the pathway of the source streamline but also the scalar values attached to it such as FA or GFA. This would permit to discover dierences related to both the geometry of the streamlines and to their microstructure.

Topological changes in ber bundles

The double dieomorphic scheme presented in Chapter 4 can not take into account all the topological variations of the white matter. For instance, it can not put into correspondence circuits with intersecting or shrinking bers. In Fig. 5.1 we show on the left a template with parallel streamlines connecting homogeneously

Limitations and Perspectives

two toy grey matter structures and on the right a subject complex with intersecting streamlines connecting only specic grey matter areas. The only way to attain a correct matching would be to use a dieomorphism for each streamline. Instead, in This produces a data-term which depends entirely on the norm of the template and therefore the only way to decrease it is shrinking the template bers.

The rst problem would happen mainly when considering several ber bundles at the same time, like all the bundles of the CSTC circuits, instead the second one when there are important dierences in structural connectivity among the population under study or between the chosen initial template and the subjects bundles.

We envisage here two possible solutions to solve these problems. The rst one would be to use a dierent deformation framework: the one of metamorphosis (Trouvé and Younes, 2005;[START_REF] Holm | The Euler-Poincaré theory of metamorphosis[END_REF][START_REF] Younes | Computing metamorphoses between discrete measures[END_REF].

It combines a dieomorphic mapping with a topological change of the template.

This would permit to make disappear the shrinking bers from the template or to change the arrangement of its streamlines in order to take into account the intersecting streamlines. Moreover, in this framework it is possible to compute population statistics based on the initial momenta parametrising the deformations, as shown in [START_REF] Richardson | Metamorphosis of images in reproducing kernel Hilbert spaces[END_REF].

A second (and more radical) solution would be to discard the pathway of the streamlines keeping only the information about the structural connectivity. Given two ROIs (Region Of Interest), which could be the entire putamen and the entire left hemisphere of the cortex as in Chapter 4, we would compute a probability map at each vertex of the rst ROI indicating the probability to be linked with any vertex of the second ROI. This probability map would describe the structural connectivity.

Then we could use functional currents (or functional varifolds) (Charon and Trouvé, 2013a) to model the grey matter surfaces where the functional signal would be the probability map. Eventually, one could use a kind of double deformation to modify rst the probability maps and then both the shape of the structures and their functional signal. This would avoid all topological problems about ber bundles but it would discard the information about the pathway of the streamlines, which might be important, and it would depend on the choice of the ROIs.

Clinical variables

A last perspective is about the integration of clinical variables in the proposed framework. A clinical variable is a scalar value which quanties a specic characteristic of a pathology aecting a patient like the severity of the symptoms. An example is the scalar YTGTSS/50 which quanties tic severity.

Clinical variables could be used in dierent ways. We could estimate the PCA modes that best explain the clinical variables with linear regression analysis and model reduction methods such as Akaike information criterion (AIC). Furthermore,

Limitations and Perspectives

template and the clinical variables, which means understanding how the template changes as the clinical variable varies. A similar approach was proposed for the quantication and prediction of cardiac remodelling in [START_REF] Mansi | A Statistical Model for Quantication and Prediction of Cardiac Remodelling: Application to Tetralogy of Fallot[END_REF].

Eventually, we could also integrate shape and clinical information into a single framework by using manifold learning techniques [START_REF] Aljabar | A Combined Manifold Learning Analysis of Shape and Appearance to Characterize Neonatal Brain Development[END_REF]. Pairwise distances between couples of subjects could be dened by the amount of deformation needed to warp the rst one onto the second one, i.e. the norm of the initial velocity eld parametrising the deformation. In our case, we could use the two dieomorphisms together or separately. These distances could then be used as input

1

 1 Atypical brain conguration in neurodevelopmental disorders . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Gilles de la Tourette syndrome . . . . . . . . . . . . . . . . . 3 1.3 Driving hypothesis about GTS pathophysiology . . . . . . . 3 1.4 The role of structural and diusion imaging . . . . . . . . . 4 1.5 A holistic approach . . . . . . . . . . . . . . . . . . . . . . . . 6 1.6 Mesh-based morphometry . . . . . . . . . . . . . . . . . . . . 8 1.6.1 Computational models for complexes of 3D meshes . . . . . . 9 1.6.2 Statistical shape analysis of complexes of 3D meshes . . . . . 10 1.6.3 Structural connectivity analysis . . . . . . . . . . . . . . . . . 13 1.7 Contributions and Manuscript Overview . . . . . . . . . . . 14 1.8 List of Publications . . . . . . . . . . . . . . . . . . . . . . . . 18 1.1 Atypical brain conguration in neurodevelopmental disorders Neurodevelopmental disorders such as autism (Belger et al., 2011; Tye and Bolton, 2013), attention decit-hyperactivity disorder (ADHD) (Konrad et al., 2010; Xia et al., 2012) or Gilles de la Tourette syndrome (GTS) (Liu et al., 2013; Worbe et al., 2015) are characterised by an atypical brain development which is thought

Figure 1 . 1 :

 11 Figure 1.1: Schematic diagram illustrating the three main cortico-striato-thalamocortical circuits. Figure taken from (Krack et al., 2010).

Figure 1 . 2 :

 12 Figure 1.2: Two views of the same anatomical complex of 3D meshes representing the caudate, putamen and thalamus of the left hemisphere of the cortex in red, yellow and blue respectively. Streamlines connecting these three sub-cortical nuclei to the left hemisphere of the cortex are also shown with darker colours. Fibers stem from a deterministic tractography and they have been uniformly downsampled for visualisation purpose. B=back of the head, F= front of the head.

  ; Marsland and McLachlan (2007); Avants et al. (2008); Ashburner and Friston (2011); Vialard et al. (2011).

  Ashburner (2007); Vercauteren et al. (2008); Modat et al. (2012); Lombaert et al. (2012); Lorenzi et al. (2013). Other dieomorphic registration algorithms exist such as the one in Christensen et al. (1996) which models the template as a viscous uid whose deformations are based on the Navier-Stokes equations (non-linear partial dierential equations) or the one in Rueckert et al. (2006) based on Free-Form deformations. For a more exhaustive description about dieomorphic registration algorithms the user is referred to Sotiras et al. (2013).

2. 1

 1 . Introduction(2014) where the authors have shown the importance of the information about the relative position between separated structures to discriminate controls from patients using sub-cortical nuclei of the brain.

  Durrleman et al. (2011b); Allassonnière et al. (2007); Ma et al. (2008):

Glaunès ( 2005 )

 2005 and it can be used for both surfaces and curves. It does not require point-correspondence between structures or curve-correspondence between ensemble of tracts. This means that the same structure may have a dierent number of vertices among subjects. It is robust to mesh imperfections such as holes, spikes or tract interruptions. Moreover, it does not require a consistent orientation of the normal or tangent vectors among the population since the distance is invariant to a change of orientation of some normals or tangent vectors. Another important characteristic is that varifolds prevent what is known as the currents cancelling eect. This happens when two surface/curve elements with opposite orientation cancel reciprocally their eect on the estimation of the template. For a more detailed discussion about varifolds the user is referred to Durrleman et al. (2014); Charon and Trouvé (2013b).

  ; Durrleman et al. (2014).

  pq where the scalars c uk and c pq refer to the values obtained at the end of the projection of the varifold in {x u , ← → β k } and {x p ,

( 2 .

 2 20)ODEs are integrated with the Euler's method using 10 steps. More details can be found in the Appendix of[START_REF] Durrleman | Morphometry of anatomical shape complexes with dense deformations and sparse parameters[END_REF].

Figure 2 . 1 :

 21 Figure 2.1: Template initialisation for three brain ber bundles: left caudate, left putamen and left thalamus bundle. The top row contains the three raw initial templates where the bers of every subject present in the population have been gathered together (see 2.3.1). The second row presents the initial templates obtained at the end of the proposed initialisation process cut at the intersection with their respective sub-cortical initial templates, transparent and coloured in red.

2. 3 .

 3 Experiments sphere of the cortical surface. The three sub-cortical structures are segmented with FSL Patenaude et al. (2011) from T1-weighted images (3T). Since we want to analyse as a single structure nucleus accumbens and caudate, we merge these two segmentations together. We use the marching cubes algorithm to create the 3D meshes for all structures (AimsMesh function of Brainvisa 4.3.0). Fiber bundles come from a probabilistic tractography method using 8 seeds Perrin et al. (2005) applied to diusion-weighted images (50 directions, B-factor=1000) using an analytical Q-Ball model to estimate the local underlying orientation distribution function (ODF) Descoteaux et al. (2007). Every subject bundle is downsampled by randomly selecting 10% of its bers and it is then approximated as proposed in Durrleman et al. (2009) using a Matching Pursuit Algorithm. For more details about the acquisition or the tractography, the reader is referred to Worbe et al. (2015).

  [10 -2 -10 3 ], the number of possible combinations of the values of the parameters for the Fixed technique is 36 whereas for the hyperparameters of the Bayesian framework is 7776. This makes infeasible an extensive and complete analysis of the robustness of the results. Nevertheless, we have decided to compute 30 atlases (15 Fixed + 15 Bayesian) picking randomly the values of the hyperparameters in the range [10 -2 -10 2 ]. The used hyperparameter and σ 2 j values are listed in Table

Figure 2 . 2 :

 22 Figure 2.2: Comparison of the robustness with respect to the hyper-parameters values using the proposed algorithm (Bayesian) and to σ 2 j without automatic estimates (Fixed). In all atlas estimations based on the Bayesian method, it has been changed only one of the hyperparameters xing the others to 0.01. Surfaces are the left caudates and Bundles are the bers connecting them to the left hemisphere of the cortical surface. Figures on the left represent the norm of the dierence between the template obtained at the end of the atlas constructions and a reference template, which is the one obtained using the Fixed method with σ 2 j =0.01. Figures on the right refer to the norm of the residuals obtained at the end of the atlas procedures.

Figure 2 . 3 :

 23 Figure 2.3: Comparison of the robustness with respect to the hyper-parameters values using the proposed algorithm (Bayesian) and to σ 2 j without automatic estimates (Fixed). The template complex in every atlas construction is composed by both the left caudate and the left caudate bundle. The values of the hyperparameters for the Bayesian estimation and the values for the two σ 2 j with the Fixed method are listed in Table 2.2. Figures on the left represent the norm of the dierence between the template obtained at the end of the atlas constructions and a reference template, which is the one obtained in the 9-th case using the Fixed method. Figures on the right refer to the norm of the residuals obtained at the end of the atlas procedures.

Figure 2 . 4 :

 24 Figure 2.4: Robustness with respect to the number of points Λ of the varifold grid. Nine dierent atlases are computed using left-caudate surfaces changing only Λ within the range [1-10 5 ]. Top-left gure represents the norm of the dierence between a reference template (Λ=1) and the templates obtained at the end of the atlas constructions. Top-right gure shows the norm of the residuals during the atlas constructions. Bottom gure is about the evolution of log(|Γ α |) during the atlas procedures. The maximum number of iterations in every atlas procedure is xed to 150. Colors refer to the dierent number of grid points.

Figure 2 . 5 :

 25 Figure 2.5: Atlas construction process. From left to right, we present rst the initial template complex and the initial set of control points. Then, we show the nal template obtained at the end of the atlas construction. The top and bottom arrows point to the nal template deformed accordingly to the averages of initial momenta of controls µ (top) and patients -µ (bottom). The averages of initial momenta (µ and -µ) are shown respectively above and below the two arrows.

Figure 2 . 6 :

 26 Figure 2.6: Shape dissimilarities about the sub-cortical structures between the two groups. The colors refer to the absolute value of the dierence between the displacements from the nal template to the average congurations of patients and controls. The two average congurations are obtained by deforming the nal template along the directions given by µ and -µ. The four frames represent the same three structures from dierent points of view. Letters {c, p, t} refer to caudate, putamen and thalamus respectively.

Figure 2 . 7 :

 27 Figure 2.7: Shape dissimilarities about the ber bundles between the two groups. The colors refer to the absolute value of the dierence between the displacements from the nal template to the average congurations of patients and controls. The two average congurations are obtained by deforming the nal template with µ and -µ.

Figure 2 .

 2 Figure 2.8: PCA of the two covariance matrices of the initial momenta based on the assumption Γ c α = Γ p α . The magenta line represents the percentage of variability explained by the rst mode. The red line shows the number of components needed to explain 90% of the total variability. The green line indicates the magnitude of the eigenvalues and the blue line shows the explained variability using all the previous components.

Figure 2 . 9 :

 29 Figure 2.9: Second mode of variation at ± 1.5 standard deviation (σ) of two PCAs based on the two estimated covariance matrices of the initial momenta Γ c α and Γ p α . Colors refer to the displacement from the nal common template shown in the middle of the gure in gray. Arrows indicate the areas of the template complex which vary in a dierent way between the two groups.

  ous works Durrleman et al. (2014); Avants and Gee (2004); Ma et al. (2010), the authors proposed to use the sample covariance matrix computed at the end of the atlas procedure. Since the number of subjects is usually smaller than the number of deformation parameters, the sample covariance matrix might be ill-conditioned.

Contents 3 . 1 54 3. 2 . 1 55 3. 2 . 2 55 3. 3

 3154215522553 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54 3.2 Related Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . Compact representation of streamlines . . . . . . . . . . . . . Compact representation of the whole bundle . . . . . . . . . . Our contribution . . . . . . . . . . . . . . . . . . . . . . . . . . 56 3.4 Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 3.4.1 Weighted currents . . . . . . . . . . . . . . . . . . . . . . . . 57 3.4.2 Approximation scheme . . . . . . . . . . . . . . . . . . . . . . 59 3.5 Experiments and Results . . . . . . . . . . . . . . . . . . . . . 64 3.5.1 Materials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 3.5.2 Numerical aspects . . . . . . . . . . . . . . . . . . . . . . . . 64 3.5.3 Weighted prototypes representation . . . . . . . . . . . . . . . 66 3.5.4 Qualitative evaluation of the parameters inuence . . . . . . 67 3.5.5 Registration-based evaluation of the algorithm . . . . . . . . 68 3.6 Discussion and Conclusions . . . . . . . . . . . . . . . . . . . 70 3.A Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71 3.A.1 Interpretation of τ . . . . . . . . . . . . . . . . . . . . . . . . 71 3.A.2 Modularity based on weighted currents . . . . . . . . . . . . . 71 Chapter 3. Parsimonious approximation for white matter ber bundles based on weighted currents 3.1 Introduction Tractography Conturo et al. (1999); Mori and van Zijl (2002) from diusion-weighted magnetic resonance imaging (DW-MRI) Basser et al. (1994) is the only non-invasive technique capable to trace in vivo the wiring architecture of the human brain white matter. It is widely employed for both clinical (i.e. stroke Mukherjee (2005), surgical procedures Ciccarelli et al. (2008)) and research purposes (i.e. Alzheimer's disease, schizophrenia Kubicki et al. (2007)). The 3D polylines stemming from tractography algorithms, called streamlines or bers, are only estimates of the trajectories of large groups of neural axons. Streamlines are traced from points inside a starting voxel, called seeds, and they are constituted of segments connecting neighbouring voxels. The direction of these segments is dened by a local diusion model (i.e. tensor, Q-ball Fillard et al. (2011)) computed at each voxel and by a tractography method: deterministic or probabilistic Fillard et al. (2011). Deterministic algorithms produce segments which follow the principal direction of the local diusion model whereas probabilistic ones use randomly perturbed versions of the main direction. Seeds are usually placed in every voxel of the white matter (whole-brain tractography) and the resulting streamlines can then be divided into dierent ber bundles based on clustering algorithms or starting/ending Regions of Interest (ROI). Fiber bundles may be then decomposed into fascicles which are groups of bers with a similar pathway and whose extremities are close to each other, connecting therefore the same functional territories. Fiber bundles are dicult to analyse both qualitatively and quantitatively due to their considerable number of streamlines. The size of a bundle can make intractable processes such as clustering O'Donnell and Westin (2007), registration Siless et al. (2012), atlas construction Durrleman et al. (2011a) or shape analysis Corouge et al. (2004) due to the great computational load and execution time. It can also complicate the rendering and the visualisation, thus limiting possible clinical applications.

  consists of dening point-to-point correspondence among streamlines parametrising them as sets of points O'Donnell et al. (2012) or with cubic B-splines Corouge et al. (2006) for instance. This technique eases the computations but it can be applied only if streamlines have a similar length and the denition of corresponding points can be very challenging. Other authors proposed to characterise a streamline using only its extremities Brun et al. (2003) or its connectivity signature Tunç et al. (2013), namely the probabilities to be linked to a dened set of ROIs. These methods have been used for clustering, visualisation and interpretation purposes but they do not take into consideration the shape of the streamlines which is important for registration and morphometry. Conversely, dierent authors proposed to evaluate only the local pathway of the bers without taking into consideration their extremities. A rst example is the group of methods based on Fourier descriptors Batchelor et al. (2006); Chung et al. (2010), which result in a concise parametrisation useful for clustering and shape analysis. In these models the number of descriptors needs to be xed though and the optimal number depends on the length and shape of the streamlines which might vary even within a single bundle. Lately, other authors proposed to represent a streamline as a blurred indicator function modelled as a Gaussian process Wassermann et al. (2010); Liu et al. (2012). This representation permits

  bundle. The most common strategy is to employ a centroid which represents the mean trajectory of the bers of the bundle O'Donnell and Westin (2007); Corouge et al. (2006); Guevara et al. (2011); Garyfallidis et al. (2012). Other authors have also employed isosurfaces to represent the spatial variation of the centroid within the bundle Maddah et al. (2007). This representation can be used only for tubularshaped bundles that can be modelled as convex envelopes. Other bundles, such as Chapter 3. Parsimonious approximation for white matter ber bundles based on weighted currents the corpus callosum and the rostral part of the corticospinal, have a dierent topology and they are dened as sheet-like bundles. In Maddah et al. (2011) the authors proposed to represent those bundles as 3D surface meshes whereas in Yushkevich et al. (2008) it was suggested to use deformable medial models (cm-reps). In both cases, the medial surface representations are employed only for visualisation and clustering and to provide statistics about diusion coecients. A dierent representation, which can be employed for any kind of bundle, is the tract probability map Hua et al. (2008); Bürgel et al. (2006); Wassermann et al. (2010). It indicates the probability of a voxel to belong to a given bundle. This method is very concise but it is not based on a geometrical primitive and it has been used for visualisation, interpretation and clustering. A last example is the sparse representation based on the matching pursuit algorithm for currents presented in Durrleman et al. (2011a). In the framework of currents Vaillant and Glaunès (2005) a bundle is considered as a single mathematical object composed of disconnected oriented points which model the local orientation of the streamlines. The approximation presented in Durrleman et al. (2011a) represents a bundle with a sparse set of oriented points. This representation is very concise but it has the drawback to approximate accuratelyonly the areas of the bundle characterised by a high density of streamlines, like the central mass of the bundle. Thus, the small fascicles may not be well approximated.

currents

  [START_REF] Vaillant | Surface Matching via Currents[END_REF] the fact that it does not need neither pointto-point nor streamline-to-streamline correspondences. Nevertheless, it requires to dene the starting and ending point of each streamline. The uncertainty of tractography algorithms near the grey/white matter boundaries is taken into account by a Gaussian smoothing kernel. Prototypes are visualised as tubes whose constant radii are proportional to their weights. This concise representation can be easily combined in multi-object studies with other geometric primitives such as surface 3.4. Method meshes modelled as landmarks, currents or varifolds Gori et al. (2013a); Durrleman et al. (2014).

Figure 3 . 1 :

 31 Figure 3.1: Two 2D streamlines X and Y are compared using both the framework of currents and the one of weighted currents. In the three gures only the position of the extremities change, the overall pathway remains almost unchanged. The framework of weighted currents is more sensitive to the distance between the extremities which explains why the two streamlines X and Y are almost orthogonal in the last gure on the right. The bandwidths of all kernels employed in both frameworks are equal to 5.

Figure 3 . 2 :

 32 Figure 3.2: Streamlines that have an angle smaller than 45 degrees with the red one using currents (green, #118) and weighted currents (blue, #8). Green streamlines are more spread than the blue ones, connecting anatomical locations far from the ones of the red ber. The concept of similarity in the framework of weighted currents is more stringent than using usual currents.

2 W 1 j∈F 2 Figure 3 . 3 :

 21233 Figure 3.3: Example of weighted prototype visualised as a red tube. It approximates the fascicle of streamlines coloured in green. The streamline chosen as prototype is the central axis of the tube and it is coloured in black. The radius of the tube is proportional to the weight of the prototype.

Figure 3 . 4 :

 34 Figure 3.4: Visual explanation of the bundle subdivision into fascicles (modes).Dots and crosses represent the streamlines of a ber bundle in the space of weighted currents. The green cross is the most similar streamline to the average of the bundle.It is far from almost all the other streamlines and it could be considered as an outlier.The three circles represent the fascicles (modes) composed by similar bers. The red crosses are the prototypes of the fascicles considered independently. These bers are more representative than the green cross and they better approximate the ber bundle.

  Diusion weighted scans are acquired with sequences of 50 directions with a B-factor of 1000 and a voxel size of 2x2x2 mm 3 . We use the Spherical Deconvolution Transform (SDT) model Descoteaux et al. (2009) to estimate the local underlying orientation distribution function (ODF). Whole brain connectivity is then inferred within an anatomy-based tractography mask Guevara et al. (2011) using both a deterministic (1 seed per voxel) and a probabilistic (8 seeds per voxel) tractography algorithm available in BrainVISA/Connectomist-2.0 Perrin et al. (2005). In this chapter, we consider three distinct ber bundles connecting the left hemisphere of the cortical surface to the left thalamus, putamen and caudate respectively. We extract them from both the deterministic and probabilistic whole brain tractography as explained in Worbe et al. (2015). All bundles include also the commisural bers which are truncated at the interhemispheric ssure. All the other bers are cut at the intersection with their respective sub-cortical structure and at the border between white and gray matter of the cortex. Sub-cortical structures are segmented with FSL Patenaude et al. (2011) from 3D T 1 -weighted images (voxel size: 1x1x1 mm 3 ) and we merge the segmentations of nucleus accumbens and caudate in order to consider them as a single structure. The 3D meshes are created using the marching cubes algorithm of BrainVISA v4.4.0. The cortical surface is segmented using FreeSurfer v5.3 Fischl et al. (2004). More information about the acquisition and the preprocessing of both T1-weighted and diusion-weighted images can be found in Worbe et al. (2015).

Figure 3 . 5 :

 35 Figure 3.5: Weighted prototype approximations of two probabilistic bundles: a cortico-putamen and a cortico-thalamus. As it is possible to notice, our approximation does not alter neither the global shape of the bundle nor the densities of the endpoints onto the cortical surface. We use: γ=0.13, λ g =7mm, λ a =5mm and λ b =10mm.

Fig. 3 .

 3 Fig.3.5 and in Fig.3.6 composed of 80.000 and 35.674 were of 150 and 19 minutes respectively.

Figure 3 . 7 :

 37 Figure 3.7: Results of the fascicles detection step applied to the bundle shown in Fig.3.6 using the two sets of kernels' bandwidths. The number of fascicles are respectively: 65 (left) and 35 (right). Colours are chosen randomly.

Fig. 3 . 7 .

 37 Fig.3.7. The fascicles of the rst set of parameters are smaller in size and greater in number with respect to the fascicles based on the second row of parameters.

Figure 3 . 8 :

 38 Figure 3.8: On the left: registration between the green source bundle and the red target bundle approximated with weighted prototypes. On the right: original ber bundles and the transformations based on the approximation at γ=0.13 and with the sub-sample of 5000 streamlines. Black arrows highlight the areas where the alignment is more noticeable.

  λ a and λ b ) are considered as parameters xed by the user. Their values are chosen by looking at how much the streamlines fan out when approaching to the boundary between white and grey matter. Streamlines deviate more when they are close to the cortex than to the sub-cortical nuclei, which explains why λ b is always greater than λ a in our experiments. It would be of interest to automatically estimate these parameters taking into consideration the type of bundle, the SNR of the diusion image, the tractography algorithm and the diusion model. Another parameter xed at the beginning of the algorithm is the approximation level γ. It denes the stopping criteria and the value (1-γ) is the minimal percentage of the norm of the fascicle explained by the prototypes. In fact, thanks to the triangle inequality, we can rewrite the stopping criteria obtaining

2 N i=1 S i 2 W

 22 sum of the weights of all the edges joining only the vertices of module c, S c = i∈c N j=1 S i , S j W * is the sum of the weights of the edges between the vertices in c and all the N vertices in the graph and m=1/* is the sum of the weights of all edges in the graph.Substituting these equations in Eq.3.4 and noting that m is a constant term and that it can be rewritten as 2m= i∈c S i , S j W * , one obtains Eq.3.3. Chapter 4. A deformation framework to unify morphometry and structural connectivity analysis 4.1 Introduction The pathophysiology of neurodevelopmental disorders such as autism, attention decit-hyperactivity disorder (ADHD) and Gilles de la Tourette syndrome (GTS) involves morphological alterations of the cortico-basal ganglia and cortico-thalamus neural circuits (Tye and Bolton, 2013; Konrad et al., 2010; Worbe et al., 2015).These networks are composed of neural projections linking particular areas of the cortical surface to specic parts of the sub-cortical nuclei. Abnormalities can aect: i) the shape of the single components of the circuits, including both gray matter structures and white matter tracts, ii) the relative position among gray matter structures and iii) the structural connectivity, namely the areas where gray matter integrate white matter tracts. Most of the studies present in the literature focus only on the rst or last point. Few of them propose to analyse the rst two points together and, to the best of our knowledge, no method has been proposed to tackle the three points in a unied framework.

Chapter 4 .

 4 A deformation framework to unify morphometry and structural connectivity analysis cost function and using a gradient descent scheme. The data-term depends only on the deformed template complex obtained at the end of the second dieomorphism. Looking again at the example in Fig.4.1, the rst dieomorphism would make the ber bundle slide from the left gyrus to the right gyrus. The second dieomorphism would then modify the shape of all structures producing an accurate matching. The rst dieomorphism would capture therefore the changes in structural connectivity whereas the second one would recover the global morphological dierences.

Figure 4 . 1 :

 41 Figure 4.1:A template and a subject shape complex composed by a pseudo cortex, divided into a black and green gyrus, a blue pseudo sub-cortical nucleus and a red pseudo ber bundle. A single dieomorphism could not put into correspondence all the structures and capture the dierences in structural connectivity. The points within the violet circle in the template complex would be matched either to the black gyrus of the subject shape complex or to its ber bundle. A double dieomorphism would rst move the ber bundle from the left to the right gyrus and then it would change the shape of all the structures, producing an accurate matching and capturing also the variation in structural connectivity.

3 .

 3 Fiber bundles are approximated with a parsimonious representation of weighted streamlines prototypes. Prototypes are chosen among the streamlines of the ber bundle and they represent groups of similar streamlines. This approximation is controlled by the metric of weighted currents.We propose to model gray matter structures as varifolds(Charon and Trouvé, 2013b) or landmarks if correspondences across subjects are available. The framework of varifolds is the non-oriented extension of the one of currents. It does not require neither point-correspondence nor a consistent orientation of the normal vectors among the population. Moreover, it prevents the currents cancelling eect which happens when two surface cells modelled as usual currents and with opposite orientation (i.e. a bump) cancel reciprocally their eect on the estimation of the template. A structure which is typically aected by this eect is the thin tail of the caudate nucleus. The atlas is estimated within a Bayesian framework based on a generative statistical model similar to the one presented in Chapter 2 (Gori et al., 2013a; Allassonnière et al., 2007) and adapted to double dieomorphisms. Every shape is modelled as a sum between a double dieomorphism of the template complex and a noise. The noise and the deformation parameters of both dieomorphisms are modelled as Gaussian random variables centred at zero. We use priors similar to the ones in Allassonnière et al. (2007) which enable us to automatically estimate the noise variance of each structure and a well-conditioned covariance matrix of the deformation parameters for both dieomorphisms. The noise variances balance the importance of the structures in the cost function with respect to the other structures and with respect to the regularity terms of the deformations. As demonstrated in Chapter 2 Chapter 4. A deformation framework to unify morphometry and structural connectivity analysis the results are rather sensitive to their values and their automatic estimate becomes even more important increasing the number of structures.

  the triples {y l , f a m , f b n }. The matrix K Q,Λ j has a dimension of [Λ j xΛ j ].The squared norm of the projection of the weighted current dened in Eq.4.12 is:

(

  LDDMM) framework based on the control point formulation presented in Durrleman et al. (2011b). For every subject i, both φ W i and φ All i are dened as the last deformations of two ows of dieomorphisms {φ W it } t∈[0,1] and {φ All it } t∈[0,1]

  [START_REF] Gori | Bayesian Atlas Estimation for the Variability Analysis of Shape Complexes[END_REF] A sketch of the cascade of dieomorphisms is shown in Fig.4.2.

Chapter 4 .

 4 A deformation framework to unify morphometry and structural connectivity analysis

Figure 4 . 2 :

 42 Figure 4.2: Sketch of the deformation of the template complex based on the proposed cascade of two dieomorphisms φ W and φ All . We omit the subject index i for clarity purpose.

  as the ones of T and L All i respectively.

  values of θ All i relative to the grey matter structures. A sketch of the optimisation procedure to update all parameters is shown in Fig.4.3. It should be noted that the ow of information goes in the opposite sense, namely from the second global dieomorphism All to the rst white one W , with respect to the one shown in Fig.4.2 to deform the template complex T . More details about the computations can be found in Sec.4.A.2.

Figure 4 . 3 :

 43 Figure 4.3: Sketch of the optimisation procedure for the template complex and the deformation parameters. We omit the subject index i for clarity purpose.

Chapter 4 .

 4 A deformation framework to unify morphometry and structural connectivity analysis4.2.6 Initialisation of the parametersThe use of a gradient descent scheme to optimise the atlas parameters requires an initialisation. Control points of both dieomorphisms are initialised as a regular lattice covering the entire ambient space with a inter-points distance equal to the bandwidth of the dieomorphic kernel K D . Momenta are initialised to zero, namely no deformations. The third parameter is the template complex which is composed of two kinds of meshes: 3D surfaces and bundles of streamlines approximated as weighted prototypes. The initial template should not be chosen among the meshes of the subjects since it would inherit all the imperfections due to segmentation errors or data noise, biasing consequently the atlas estimation.

  3.0 which allows to have vertex-correspondence between subjects. The putamen and caudate are segmented with FSL (Patenaude et al., 2011) from the T1-weighted images. Since we want to analyse as a single structure nucleus accumbens and caudate, we merge these two segmentations together. We use the marching cubes algorithm to create the 3D mesh (AimsMesh function of Brainvisa 4.3.0). The ber bundles result from a deterministic tractography algorithm (1 seed per voxel) available in BrainVISA/Connectomist-2.0 (Perrin et al., 2005). We use the Spherical Deconvolution Transform (SDT) model (Descoteaux et al., 2009) to estimate the local underlying orientation distribution function (ODF). Whole brain connectivity is then inferred within an anatomy-based tractography mask (Guevara et al., 2011) from which the ber bundles are extracted as explained in[START_REF] Worbe | Altered structural connectivity of cortico-striatopallido-thalamic networks in Gilles de la Tourette syndrome[END_REF]. The bundles include also the commisural bers which are truncated at the interhemispheric ssure. All the other bers are cut at the intersection with their respective sub-cortical structure and at the border between white and gray matter of the cortex. For more details about the acquisition, preprocessing of both T1weighted and diusion-weighted images or the tractography, the reader is referred to[START_REF] Worbe | Altered structural connectivity of cortico-striatopallido-thalamic networks in Gilles de la Tourette syndrome[END_REF].

Fig. 4 .

 4 Fig.4.4 we show the results where V and W C indicate the number of points of the Varifold and Weighted Currents grid respectively. The rst row represents the evolution of the residuals of both structures during the atlas construction. Similarly, the second row represents the evolution of the logarithm of the determinant of the two covariance matrices, namely log(|Γ W α |) and log(|Γ All α |). The last row describes instead the point-wise Euclidean dierence || • || 2 between the estimated template and the initial one.

Figure 4 . 4 :

 44 Figure 4.4: Evaluation of the robustness of the results with respect to the number of points of the varifold (V) and weighted current (WC) grids. Colours refer to 16 atlases computed using dierent values of grid sizes. Each atlas is based on 8 anatomical complexes composed of left caudate, modelled as varifold, and left caudate bundle, modelled as weighted current. The maximum number of iterations for all atlases is 100.

4. 3 .

 3 Experiments and Results

Figure 4 . 6 :

 46 Figure 4.6: Toy example of double dieomorphic atlas construction. At the top we present the shape complexes of two dierent populations. The intra-variability of each group is smaller than the inter-variability. Subjects of A have a bigger cortex than subjects in B and both the shape of the sub-cortical nucleus and the structural connectivity are dierent. In the middle we show the initial template which is updated during the atlas construction. The nal template is presented at the bottom and it is deformed at ±σ along the rst modes of two Principal Component Analysis (PCA) computed with the covariance matrices of the momenta of the two dieomorphisms. Due to the great inter-variability between the two groups, the rst mode of the white PCA (rst dieomorphism) presents the main dierences in structural connectivity between the two populations. The rst mode of the global PCA (second dieomorphism) shows instead the global morphological variations between the two groups.

Figure 4 . 7 :

 47 Figure 4.7: Example of template-to-subject registration using both a single and a double dieomorphism. Due to the variation in structural connectivity between the template complex (blue, green and red) and the subject complex (violet), a single dieomorphism can not put into correspondence all structures. Instead, a double dieomorphism can rst make sliding the bers onto the grey matter structures, showing the variation in structural connectivity, and then it can correctly register all structures. Orange arrows indicate the changes in structural connectivity.

Figure 4 . 8 :

 48 Figure 4.8: Main variations in structural connectivity of an atlas composed of 10 controls and 10 patients. The nal template bundle, shown in the middle, is deformed at ±3 std along the rst mode of a PCA based on the initial momenta of the white dieomorphisms. Colours of the left and right columns indicate the densities of the endpoints of the template bundle onto the xed grey matter structures. Colours in the middle column highlight instead the displacement of the template bundle at the end of the deformations. Arrows point to the most variable areas.

Figure 4 . 9 :

 49 Figure 4.9: Main global morphological variations of an atlas composed of 10 controls and 10 patients. The nal templates of the grey matter, shown in grey in the middle row, are deformed at ±3 std along the rst and second mode of a PCA based on the initial momenta of the global dieomorphisms. Colours on the left highlight the displacement of the template cortex at the end of the deformations.

Figure 4 . 10 :

 410 Figure 4.10: Main variations in structural connectivity of an atlas estimated using 27 healthy controls. Each row shows the nal template bundle deformed at ±3 std along the rst three modes of a PCA based on the initial momenta of the white dieomorphisms. Colours refer to the displacement (in mm) of the template bundle at the end of the deformations.

Figure 4 . 11 :

 411 Figure 4.11: Main variations in structural connectivity of an atlas estimated using 30 GTS patients. Each row shows the nal template bundle deformed at ±3 std along the rst three modes of a PCA based on the initial momenta of the white dieomorphisms. Colours refer to the displacement (in mm) of the template bundle at the end of the deformations.

Figure 4 . 12 :

 412 Figure 4.12: Bootstrap analysis of 1000 iterations performed on the top of a LDA with a leave-one-out cross validation. Each sample of the histogram represents the classication score obtained using 44 patients chosen randomly (with replacement) among all the sub-groups and 22 xed controls. Red and green lines show the average and the 95% condence interval respectively.

  Both dieomorphisms are parametrised by a dierent set of control points. The number of control points denes the dimension of the initial momenta. These are the features used in statistical analysis such as PCA or classication. It would be of interest investigating how a reduction of the number of control points would inuence both the atlas construction and the statistical analysis. In[START_REF] Durrleman | Morphometry of anatomical shape complexes with dense deformations and sparse parameters[END_REF], the authors used a single-dieomorphic atlas construction method similar to the one we proposed. They demonstrated that the statical performance in a classication augments by decreasing the number of control points until a certain threshold. It seems therefore reasonable to expect the same behaviour for the proposed method.

Chapter 4 .

 4 A deformation framework to unify morphometry and structural connectivity analysis 4.A Appendix 4.A.1 Data-term and gradient of weighted prototypes Given two sets of weighted prototypes: N i=1 τ i F i and M p=1 w p T p where F i = L i n=1 δ α in (x in ,f i ,t i ) and T p = Zp l=1 δ β pl (y pl ,vp,up) we rst present how to compute the squared norm of their dierence || N i=1 τ i F i -M p=1 w p T p || 2

τwτ

  i τ j k a (f i , f j )k b (t i , t j ) F i , F j U * p w d k a (v p , v d )k b (u p , u d ) T p , T d U * i w p k a (f i , v p )k b (t i , u p ) F i , T p U * (4.32)where U * is the space of usual currents, namelyF i , F j U * = L i n=1 L j m=1 α T in k g (x in , x jm )α jm , T p , T d U * = Zp l=1 Z d o=1 β T pl k g (y pl , y do )β do , F i , T p U * = L i n=1 Zp l=1 α Tin k g (x in , y pl )β pl .

  Here we consider only the variations related to the data term D i [T All i1 ]. The Chapter 4. A deformation framework to unify morphometry and structural connectivity analysis variations of the regularity terms can be found in Sec.4.2.5. It results: LAll i (t) = (d L All it F All i (t)) T δL All it δL All i (0) = δL All i0

t 0 d∂ 0 ∂ 0 V

 000 L(u) F (u)du δL 0 = R 0t δL 0 (4.34)The last two ODEs are instead linear with source term and their solution is given by:T (u) G(u)du ∂ L(u) G(u)δL(u)du + exp t T (s) G(s)ds δT 0 = t ut ∂ L(u) G(u)δL(u)du + V 0t δT 0 (4.35)where, as in[START_REF] Durrleman | Morphometry of anatomical shape complexes with dense deformations and sparse parameters[END_REF], we denote: R st = exp(t s d L(u) F (u)du) and V st = exp( t s ∂ T (u) G(u)du).We omit the upper indices W and All and the index i for clarity purpose. Substituting Eq.4.34 into Eq.4.35 and noting that δT All,W (0) = δT W (1), we obtain:

δT 0 V 1 0∂ 1 0∂

 011 All,W (t) = t All ut ∂ L All (u) G All (u)R All 0u du δL All 0 + V All 0t δT W 0 + (4.36) V All 0t T W (s) G W (s)δT W (s)ds + V All 0t L W (s) G W (s)δL W (s)ds

  4 a new dieomorphic scheme based on a cascade of two dieomorphisms. This permits to model variations in structural connectivity with the rst dieomorphism and global anatomical changes with the second one. The parameters of each dieomorphism are assumed to follow a dierent and independent Gaussian distribution centred at zero. A double dieomorphic atlas construction embedded in the Bayesian framework presented in Chapter 2 results in a template complex, showing the common anatomical characteristics of the population under analysis, and two well-conditioned covariance matrices, describing the organisational and morphological variability captured by the double dieomorphisms.

Fig. 5 .

 5 Fig.5.2 we show an example of shrinking bers. This happens when the bers of the template, shown on the left, do not have correspondent bers on the subject bundle, shown on the right, within an area equal to the dieomorphic kernel bandwidth.

  of a manifold learning algorithm together with the clinical information like in Wolz et al. (2011) or Fiot et al. (2012). In the resulting low-dimensional representation of the data we could then use classication algorithms such as Support Vector Machine (SVM).

Table 2 .

 2 1: Angles in degrees between the rst 5 modes of variations of the group of controls and the rst ve of the group of patients. The average angle between all

			1	2	Patients 3	4	5
			1 25.71 82.49 87.10 87.77 83.32
		Controls	2 87.70 60.39 72.57 83.62 80.30 3 80.60 62.59 67.17 84.39 79.76 4 86.74 75.15 82.60 88.98 77.35
			5 83.70 81.24 89.51 85.39 83.80
	modes is 85.39	• .

Table 2

 2 .2 the values of the xed parameters and hyperparameters used in the experiments shown in Fig.2.3.

	Chapter 2. A Bayesian framework for multi-object atlas construction

Table 2 .

 2 2: Values of σ 2 for the Fixed method and hyperparameter values for the Bayesian algorithm used in the robustness analysis shown in Fig.2.3.This chapter has been submitted to IEEE Transactions on Medical Imaging and it has been partly published in[START_REF] Gori | A Prototype Representation to Approximate White Matter Bundles with Weighted Currents[END_REF]. The C++ code about the computational model of weighted currents will be integrated to the software deformetrica.

	Chapter 3

Table 3 .

 3

	Chapter 3. Parsimonious approximation for white matter ber bundles based on weighted currents
		Deterministic
	Bundle	Cortex-caudate	Cortex-putamen	Cortex-thalamus
	N Streamlines 17079±4881	28056±5247	28371±6806
	N Prototypes	344±58	409±55	341±87
	Compression	97.85 %	98.49 %	98,77 %
	Table 3.1: Average compression ratio (%) of all the bundles obtained from the
	deterministic tractography algorithm
		Probabilistic
	Bundle	Cortex-caudate	Cortex-putamen	Cortex-thalamus
	N Streamlines 75389±4646	78125±2223	68640±8568
	N Prototypes	1182±358	1411±393	1000±319
	Compression	98,41 %	98,19 %	98,54 %

2: Average compression ratio (%) of all the bundles obtained from the probabilistic tractography algorithm

Table 3 . 3

 33 Error Computational Time

	γ=0.05	9.79e+08	18h 35min
	γ=0.13	9.88e+08	5h 51min
	γ=0.25	1.42e+09	4h 34min
	γ=0.40	1.75e+09	2h 39min
	1000 streamlines	9.96e+08	23h 12min
	5000 streamlines	9.94e+08	547h 32min
	8000 streamlines	-	∼ 1120h
	15000 streamlines	-	∼ 4484h
	30000 streamlines

-∞ : Registration error and computational time using dierent approximation levels and samples of the original bundles.

  4.2. Methodsa Posterior (MAP) estimator where the hidden variables are actually considered as parameters to be optimised. The whole set of parameters is therefore: {T ,

	Γ All α , Γ W α , {σ 2 j } j=1,...,M , {α All i0 } i=1,...,N , {α W i0 } i=1,...,N , c All 0 , c W 0 }. We use a gradient
	descent scheme for {T , {α All i0 } i=1,...,N , {α W i0 } i=1,...,N , c All 0 , c W 0 }. Instead the other
	parameters, namely Γ All α , Γ W α , {σ 2 j } j=1,...,M , have closed form solutions thanks to the
	use of conjugate priors. They result:

  Chapter 4. A deformation framework to unify morphometry and structural connectivity analysis This brings to another question which is how to choose the position and the number of control points. A possible solution was presented in Allassonnière et al. (2015) by integrating the selection of the best control points in the minimisation process. The authors proposed to start from a regular grid which is trimmed by keeping only the control points that participate to the template-to-subject deformations of the whole population. The idea is brought from the Group LASSO Literature and it seems to work well in practice. It would be of interest to integrate it to the proposed model. In this chapter, like in Chapter 2, we approximate the conditional distribution of the E step in the optimisation procedure with a Dirac distribution at its mode. It was demonstrated in Allassonnière et al. (2007) that this algorithm might not converge when the SNR (Signal to Noise Ratio) is low. A possible solution would be to employ sampling techniques such as Markov Chain Monte Carlo (MCMC) methods (i.e. SAEM-MCMC) which would converge whatever the quality of the data. Moreover, this would also permit to estimate the entire posterior distribution of the deformations, and not only its mode, which contains important information such as the registration uncertainty (Risholm et al., 2013; Simpson et al., 2012).

Another perspective regards the initialisation of the template of ber bundles. In Chapter 2 and in this chapter we proposed two dierent methods. In both cases, we rst combined the original streamlines of the bundles of the subjects in a single raw template. Then, we approximated it either by looking for the most representative streamlines (Sec.2.2.5) or by representing it with few scattered weighted prototypes (Sec.

4.2.6

  4.A. Appendixm ks = f k ), results: k a (f i , f k )k b (t i , t k ) k a (v p , f k )k b (u p , t k )If p ks = t k one needs to add the following to the previous result:

	2τ k	N i=1	τ L i n=1	±α in exp(-	||x ks -x in || 2 λ 2 g	) -α T ks α in (	x ks -x in λ 2 g	) exp(-	||x ks -x in || 2 g λ 2	) -
	2τ k	M p=1	w Zp l=1	±β pl exp(-	||x ks -y pl || 2 λ 2 g	) -α T ks β pl (	x ks -y pl λ 2 g	) exp(-	||x ks -y pl || 2 g λ 2	)
											(4.33)

i p

  If p ks = t k , the update rule of the point p ks is thus equal to The update rule of the rst point m ks = f k is instead

	∂|| N i=1 τ i F i -M p=1 wpTp|| 2 Q * ∂p ks	+	∂|| N i=1 τ i F i -M p=1 wpTp|| 2 Q * ∂m ks+1	otherwise if p ks = t k it is sim-
	ply . equal to ∂|| N i=1 τ i F i -M p=1 wpTp|| 2 Q * ∂p ks ∂|| N i=1 τ i F i -M p=1 wpTp|| 2 Q * ∂m ks	.

a )k b (t k , u p ) F k , T p U *

  Conclusive summary . . . . . . . . . . . . . . . . . . . . . . . . 119 5.2 Limitations and Perspectives . . . . . . . . . . . . . . . . . . 121 5.2.1 Mixture of templates . . . . . . . . . . . . . . . . . . . . . . . 121 5.2.2 Cortical surface . . . . . . . . . . . . . . . . . . . . . . . . . . 121 5.2.3 Integrating functional information . . . . . . . . . . . . . . . 122 5.2.4 Topological changes in ber bundles . . . . . . . . . . . . . . 122 5.2.5 Clinical variables . . . . . . . . . . . . . . . . . . . . . . . . . 123

	Chapter 5
	Conclusions and Perspectives
	Contents
	5.1

Chapter 4

A deformation framework to unify morphometry and structural connectivity analysis This chapter will be submitted to NeuroImage and it has been partly published in Gori et al. (2015a). The C++/CUDA code of the described double dieomorphic deformation will be soon available within the software deformetrica.

Contents

Chapter 4. A deformation framework to unify morphometry and structural connectivity analysis We show sensitivity, specicity and balanced accuracy using always 22 controls and either one of the patient sub-groups (G1,G2,G3) or combinations of them.

4.A. Appendix

Using the Fubini's theorem we can rewrite the 3 rd term as:

Instead the 4 th term becomes:

Plugging them with Eq.4.35 for δT All,G (t) into δE we obtain:

]

] we could also use Partial Least Squares (PLS) regression to automatically extract the modes that are relevant to both the shape variability and the clinical variables.

Such an approach may be useful to investigate whether an increase in the severity of tics is correlated with more pronounced atypical connectivity patterns or shape alterations in grey matter structures.

Another interesting extension would be to estimate the relationship between the