Quantitative image processing for disease study

Chapter 4 will finally present work towards the use of the previously introduced models to use them on patient data and group studies. This includes interpolation 1.3. Contributions summary 5 and averaging of DCM images, an atlas-based fiber analysis framework of a patient against controls, and new frameworks for the combined use of relaxometry and diffusion for detecting new patterns in MS patient lesions, including detecting lesions enhanced by Gadolinium (Gd) without using the contrast agent.

Foreword

My research activity in the past years has been focused on several topics all linked to better understanding the brain architecture and neurological diseases, and ultimately help with patient care. To keep it readable, this manuscript however had to leave out part of my research that has been carried out since the beginning of my career, particularly long running topics such as atlas construction, image segmentation and segmentation validation. These topics lead to great collaborations and several publications particularly with the Asclepios team (Liliane Ramus and Grégoire Malandain on multi-atlas segmentation) and the Computational Radiology Laboratory1 at Children's Hospital Boston, emphasized by an Inria associate team from 2011 to 2017. In addition, these topics lead to new methods for pediatric longitudinal brain analysis [START_REF] Legouhy | [END_REF], for Multiple Sclerosis (MS) lesions segmentation [START_REF] Karpate | [END_REF], Galassi et al. 2018] and the recent organization of a challenge workshop at the MICCAI 2016 conference on MS lesions segmentation [START_REF] Commowick | [END_REF]. I chose to focus this document on the research we (I and all PhD students, post-docs, interns and researchers I had the chance to work with) conducted on quantitative medical imaging to go towards a better understanding of neurodegenerative diseases and patient care adaptation and follow-up.

Quantitative images for disease understanding

In the field of disease diagnosis, Magnetic Resonance Imaging (MRI) has been playing for a long time a major role to provide a precise, yet non invasive, evaluation of the patient disease status. It is now used for many diseases and is a major tool for clinicians. Among other examples, stroke is one of the earliest conditions to have greatly benefitted from MRI [START_REF] Chalela | [END_REF], Warach et al. 1995]: with a combination of structural and very simple Diffusion Weighted (DW) scans, it enables the separation of patients among the different treatments available (medication or thrombectomy). Another example is epilepsy [START_REF] Kuzniecky | [END_REF] where MRI allows the detection of foci responsible for seizures enabling their removal by surgery. MRI is also useful for evaluating brain tumors where many studies have been performed [START_REF] Gordillo | State of the art survey on MRI brain tumor segmentation[END_REF] especially enabling the precise segmentation of tumors to remove them by surgery, treat them by radiotherapy or model their evolution. Finally, a class of diseases of interest in this manuscript is the class regrouping neurodegenerative diseases. Two emblematic ones are Alzheimer's disease and MS. Both diseases are still not fully understood and patient follow-up is key to the evaluation of the disease aggressiveness. For those diseases, MRI evaluation has become a crucial marker. In Alzheimer's disease, many studies have established a link between cortical thickness or general atrophy seen from MRI and the patient's status [Frisoni et al. 2010]. For MS, MRI has even become so crucial that it is part of the diagnosis criteria of the disease [START_REF] Thompson | [END_REF].

MRI however lacks specificity in its findings, i.e. it is well able to distinguish lesions or abnormalities but not to tell their specificities (tissue destruction level, etc.). For this reason, more research has been conducted to develop new, more specific, MRI sequences able to quantify the brain microstructure and its alteration: quantitative MRI techniques. Those sequences, although requiring the development of algorithms to extract relevant information, are very promising for better pathology characterization. I will quickly discuss as a starter the case of MS showing how these modalities have the potential to help the diagnosis and patient evaluation.

An example: multiple sclerosis and MRI

Multiple Sclerosis (MS) is an immune-mediated inflammatory disease of the central nervous system affecting more than 100,000 persons in France. It causes progressive myelin destruction and axonal loss (illustrated in Fig. 1.1.a) leading to increasing handicap for the patient, including walking and cognitive impairment. The disease course of MS is very variable between patients [Leray et al. 2010] (see Fig. 1.1.b and 1.1.c) and its exact causes remain largely unknown, advocating the development of imaging techniques for a better disease understanding. As mentioned above, this has already lead to great advances using MRI to provide clinicians with 1-clearly defined criteria for disease diagnosis [START_REF] Thompson | [END_REF], and 2-automatic image segmentation techniques [Danelakis et al. 2018] to count the number of lesions, their volume and evolution in time, all meaningful to evaluate the disease. The link between observations made on conventional MRI modalities such as FLAIR, T 1 -weighted or T 2 -weighted images and the disease status is unclear. This absence of correlation between clinical observations and MRI based observations has been denoted as the clinical-radiological paradox [Guttmann et al. 1995]. In particular, predicting from the beginning of the disease in which group a patient will be (on Fig. 1.1.c), is very difficult. Such an information would however be crucial to adapt the patient treatment, for example using stronger disease modifying drugs for patients at risk of a fast disease evolution.

The recent development of quantitative images in MRI offer great promises towards solving this problem. Their main difference compared to so-called conventional images is that they allow the quantification of some of the microstructure parameters i.e. specific, interpretable, properties of the underlying tissues. For example, Diffusion MRI (dMRI) quantifies the diffusion of water in multiple compartments (each related to a tissue type) and thus indirectly the microstructure of the white matter fiber bundles (axonal properties, fiber crossings) [START_REF] Panagiotaki | Compartment models of the diffusion MR signal in brain white matter: a taxonomy and comparison[END_REF], Filippi et al. 2001]. Relaxometry measures MR specific relaxation times (T 1 , T 2 ) [Tofts 2004] in multiple compartments and allows for the computation of a key component for MS: the proportion of myelin in each voxel [START_REF] Prasloski | [END_REF]]. In addition to these multi-compartment modalities, Magnetization Transfer Ratio (MTR) [Filippi & Agosta 2007] characterizes changes particularly linked to inflammation. Put together, these images could lead to unprecedented advances in the understanding of lesions specificities, their respective positions with respect to major fiber bundles and how they modify them, and help solving the clinical-radiological paradox. Moreover, researching on those images and how to process them could have great implications for many brain diseases.

Contributions summary

For all the previously mentioned reasons, we have worked on using two of these new modalities, relaxometry and dMRI for neurodegenerative diseases evaluation. This however causes several key problems of modeling, artifact correction, development of processing tools on those modalities. I will explore in this manuscript several of these challenges and the methods we proposed to tackle or reduce them. The manuscript will be split into three parts: 1-artifact correction and multi-compartment model estimation in diffusion imaging, 2-multi-compartment modeling and estimation from relaxometry, 3-processing such compartment images and designing frameworks for patient evaluation.

Diffusion imaging for white matter microstructure imaging

Chapter 2 will cover our recent developments on dMRI and particularly the progress to go beyond clinically used diffusion models such as the well known tensor model. New models, hereafter named Diffusion Compartment Model (DCM), consider the diffusion process inside a voxel as separated in several compartments, each representing diffusion in a specific tissue architecture. These models are very interesting for their interpretability. Their drawbacks however reside in their complex estimation. Moreover, diffusion images are corrupted by distortion that ought to be corrected before performing any computation on them. I present in this chapter advances on those two crucial points.

Getting insights on myelin degeneracy: relaxometry

Chapter 3 will be centered on the use of relaxometry images, able via the right estimation algorithms to provide Relaxometry Compartment Models (RCMs) i.e. models of the different compartments of water bound to either nothing (free water), myelin, or other cellular structures. Such information is of great interest to provide the myelin water fraction, to which dMRI is blind. I will show how, for different signal formation models, we define the estimation framework for obtaining robust estimates of the compartments weights. I will then illustrate preliminary results on clinical MS longitudinal data. This chapter explores our research around Diffusion MRI (dMRI) and diffusion modeling. This work has been conducted mainly with two PhD students I cosupervised: Aymeric Stamm and Renaud Hédouin. Several papers arose from this work, but particularly two main ones discussed in the following sections:

• Echo-Planar Imaging (EPI) distortion correction: Renaud Hédouin Bihan et al. 1986] measures, at each voxel location, the constrained local Brownian movement of water molecules. To measure this phenomenon, images are acquired with diffusion weighting in different directions, with different amplitudes hereafter respectively denoted gradient directions (or gradients) and b-values. On each of these DW images, the acquired intensities are directly depending on the amount of water diffusion along the gradient direction [Johansen-Berg & Behrens 2009], leading to lower signal in regions where the diffusion is high (as illustrated in Fig. 2.1 for the Cortico-Spinal Tract (CST) and corpus callosum). From these acquisitions, it is thus natural to infer a model, i.e. a Probability Density Function (PDF) in R 3 , describing the water diffusion in all directions at a given distance from its original position. In highly structured organs, such measures offer the great interest of inferring indirectly the internal structure of the organ. As an example, in the brain, the presence of highly oriented structures such as parallel axons in a fiber bundle constrains the water diffusion along their main orientation. As a consequence, the estimated model will be highly influenced by the presence of these fiber bundles and will indirectly describe them. Going further, the brain is not only composed of axons but also of a large variety of supporting cells (as illustrated in Fig. 2.2) and free water, each influencing the diffusion and thus the dMRI acquisition. A good model may thus be able to finely characterize this internal microstructure of the brain at each voxel and provide potential insights on their change over time, over individuals in a population or due to the activity of a disease.

For this reason, dMRI is a tool of choice for studying the brain microstructure. It has been widely used, both for clinical studies, e.g. in MS [Filippi et al. 2001, Werring et al. 2000, Rovaris et al. 2005], or to study the normal structure of the brain [Counsell et al. 2014, Scholz et al. 2014]. Among the large variety of challenges still at stake in this field, we will present in this chapter two main problems that were especially of importance for the goal of improving disease characterization:

• correction of susceptibility induced distortions in dMRI: images acquired through this technique have to be acquired fast to be clinically tractable. Such sequences, called EPI, suffer as a result from large anatomical distortions deforming the visual aspect of the brain. These artifacts have to be corrected for to enable a better interpretability of dMRI and to fuse its information with other modalities.

• diffusion model definition and estimation: a large variety of models may be defined from the set of DW images acquired, some requiring more gradients and b-values than others (and thus more time), some allowing to get insights into specific white matter microstructure parameters. It is a major challenge to first define the "good" model, i.e. a model describing clearly the parameters studied for a disease and that can be estimated from clinical data; and then to define the optimization procedures to properly and robustly estimate this model.

In the next sections I will detail the advances we proposed to tackle these issues:

Chapter 2. White matter microstructure from diffusion imaging artifacts correction will be studied in Section 2.2 and diffusion model estimation in Section 2.4.

Artifacts correction in diffusion weighted imaging

EPI and distortion artifacts

As mentioned above, dMRI requires fast acquisitions to reach clinically acceptable acquisition times.

To this end, EPI acquisitions are commonly used [Johansen-Berg & Behrens 2009] (and also very much used for other modalities such as functional MRI [Huettel et al. 2004]). Their high velocity comes from the fact that the image is acquired within a single repetition time (single-shot) instead of multiple shots in other classical sequences (gradient echo, spin echo...).

The high velocity of EPI acquisitions however comes at the cost of a high sensitivity to B 0 magnetic field inhomogeneities. Affected regions, often located at the tissue interfaces due to magnetic susceptibility effects, are either contracted or dilated along the Phase Encoding Direction (PED) [Jezzard & Balaban 1995]. Moreover, measured tissue intensities in these regions change due to the local transformation. Therefore the brain anatomy in EPI does not match with structural images that are much less sensitive to distortions. Such a correspondence is however necessary as a joint analysis is often performed: structural images are used to define regions of interest for fiber tracking or to extract lesions that are to be linked to brain microstructure properties. It is therefore necessary to perform EPI distortion correction, otherwise non linear anatomy mismatch between the modalities will lead to biased results.

There is therefore a growing field of approaches to solve this problem. First, early algorithms have considered the acquisition of a B 0 field map [Jezzard & Balaban 1995, Reber et al. 1998]. This map is in turn used to infer the local contractions and dilations, and to correct EPI. Other techniques have considered new sequences using point spread functions to obtain acquisitions a priori free of distortions [Robson et al. 1997, Chung et al. 2011, Zaitsev et al. 2004].

A very promising approach towards distortion correction considers the acquisition of two EPI sequences with opposite PED (e.g. one anterior-posterior and one posterior-anterior). Thanks to these additional acquisitions and through careful modeling of the distortion, images can be corrected. Moreover, it can be approximated [Vovk et al. 2007] that deformations due to distortion in successive EPI volumes are the same. Therefore, a complete dMRI volume can be corrected at the cost of only one supplementary b 0 acquisition with opposite PED. For this reason, this class of techniques has gained a lot of attention. [Voss et al. 2006] introduced an algorithm to estimate, from the two images, the correction displacement field based on cumulative intensity distributions along each line in the PED. This simple method strongly reduces the distortion, however it is sensitive to noise. The computed transformation also needs to be smoothed, leading to a trade-off between regularity and precision. [Andersson et al. 2003] used a pair of reversed EPI in con-junction with a discrete model of image formation for EPI. A registration-based method has also been proposed by [Irfanoglu et al. 2015] requiring a non distorted image such as a T 2 -weighted image (in addition to the reversed PED image) which is used as the central point where the two images with reversed PEDs are transformed.

Given the promise of this last class of techniques, and given our prior experience on registration algorithms using block-matching (see Appendix B), we have proposed with Renaud Hédouin [Hédouin et al. 2017] a new method for block-matching based distortion correction in EPI. Compared to previous approaches, we wished to account for the distortion model as early as possible through the introduction of a priori on the transformations being optimized rather than after matching through regularization (as it is done for example in [Irfanoglu et al. 2015]).

Distortion model

We assume that two images have been acquired: I F is the EPI forward image acquired with a classical PED (anterior-posterior for example), and I B is the EPI backward image acquired with a reversed PED (posterior-anterior in this case). The goal of EPI distortion correction is to estimate a distortion transformation field used in turn to recover a corrected image C. This field can also be used to correct an entire series of EPI acquired with anterior-posterior or posterioranterior PED. [Jezzard & Balaban 1995] have demonstrated that deformations due to B 0 field inhomogeneities appear mainly along the PED and are negligible in other directions. More precisely, we follow the distortion model as expressed in [Voss et al. 2006, Morgan et al. 2004] which assumes that I F and I B are generated from the theoretical true image C using a displacement field parallel to the PED:

C(x) = J T + (x)I F (T + (x)) C(x) = J T -(x)I B (T -(x)) (2.1)
where T + (x) = x + U (x) and T -(x) = x -U (x). J T + and J T -denote the Jacobian determinants of the local deformations which account for intensity changes in the distorted areas. It will lead to an increased intensity in the contracted areas and a decreased intensity in the dilated areas. U corresponds to the distortion displacement field which is parallel to the PED, e.g. if the PED is along the y-axis then

U (x) = [0 U y (x) 0] T .
It is assumed that T + and T -are opposite symmetric, i.e. that they share the same U up to a minus sign.

A block-matching strategy for distortion correction

The corrected image C or a surrogate of it is generally unknown. We therefore chose a registration approach that does not rely on it. A registration method has been introduced by [START_REF] Avants | [END_REF] allowing the deformation of two images towards their barycenter without having it directly appear in the algorithm. The idea, instead of looking for the transformation T between two images, is to seek the "half-way transformation" T 1/2 so that the two resampled images match as much as possible:

I F • T 1/2 ≈ I B • T -1/2 ≈ C (2.2)
We have adapted this idea to our distortion model and to the Block-Matching (BM) algorithm for its ability to simply and effectively incorporate constraints on the deformation field to fit the distortion model. To do so, we have first extended the general asymmetric formulation of BM to general symmetric registration (see Appendix B.4.3 for more details). Then, we have extended this symmetric framework to the distortion model by constraining the transformation to be along the PED. As a summary, we start from an initial transformation U 0 which can be null or coming from another coarse correction algorithm. Then we proceed as in Appendix B.4.3 using I B and I F as the images to register, with several modifications. At each step, we first resample the original images with the current transformation using the Jacobian determinant to modify intensities after interpolation. Then we estimate pairings between the images in the forward (A + = { Â+,1 , . . . , Â+,N }) and backward directions (A -= { Â-,1 , . . . , Â-,N }) using a BM algorithm. This matching incorporates adapted transformations, and adapted indexes of matching plausibility used in extrapolation (see Appendix B.3.2), to constrain the transformation a priori to match the distortion model. Finally, the end of each iteration incorporates a simple modification to ensure that the updated transformations at each step remain opposite symmetric by averaging the obtained deformation fields. We detail in the following the two major modifications: adapted linear transformations between blocks and adapted matching plausibility weights.

Adapted linear transformations between blocks

In other applications, the transformations A .,i sought between blocks are often 3dimensional translations (as detailed in Appendix B.2). In the case of EPI distortion, those transformations can be adapted to match a priori the expected features of the distortion at the block level and thus obtain a more robust transformation estimation. First the model assumes that distortions appear uniquely along the PED: a one-dimensional translation along the PED (modeled by one parameter t .,i ) is therefore sufficient. At the scale of the block, a single translation is however not enough to account for local contractions and dilations due to the distortion. We therefore added three parameters to the transformation definition: one for the change of scale due to the global contraction or dilation inside the block (s .,i ); and two skew components (k .,i and m .,i ) for the two directions complementary to the PED. Assuming the PED is the y-axis, A .,i is expressed as a 4 × 4 matrix:

A .,i =     1 0 0 0 k .,i s .,i m .,i t .,i 0 0 1 0 0 0 0 1     (2.3)
Note that having the PED on another axis will result in the line of parameters being displaced on the first or third line of the matrix. We have further studied in [START_REF] Hédouin | [END_REF]] the properties of these transformations which showed that there exists an analytical expression for their matrix logarithm. Therefore to speed up the algorithm, the BM search and transformation extrapolation is done in the Stationary Velocity Field (SVF) space. The BM step then amounts to estimating the four log-parameters of each block transformation to compute the set of optimal transformations Â+,i and Â-,i optimizing a similarity measure S between I F,l-1 and I B,l-1 (the two input images resampled by the current transformation at previous iteration l -1) using the BOBYQA algorithm [Powell 2009]:

       Â+,i = arg max Ã+,i S J exp( Ã+,i ) I F,l-1 • exp( Ã+,i ), I B,l-1 Â-,i = arg max Ã-,i S I F,l-1 , J exp( Ã-,i ) I B,l-1 • exp( Ã-,i ) (2.4)

Matching plausibility weights

From this set of optimal local transformations Â.,i , we then proceed to compute asymmetric transformation updates δS + and δS -using M-smoothing extrapolation as in Appendix B.3.2.2, further used to compute the update transformation. This extrapolation requires plausibility weights for the matches w .,i that provide an estimate of the confidence in the block match. In [Hédouin et al. 2017], we refined these weights to account for the uncertainty in matching along a specific direction (the PED). To do so, we use a geometric mean of two different terms. The first one is a function of the similarity at the optimal position Ŝ.,i (so that it belongs to the range [0, 1]). The second one gives an index of the local structure of the reference block along the PED. If the block structure is parallel to the PED, all tested transformations A .,i for that block would get roughly the same similarity score, introducing uncertainty in the matches. We avoid such random solutions with an index w d,i , a function of the local structure tensor inside the reference block to give a low weight to uncertain blocks and their corresponding local transformations (see [Hédouin et al. 2017] for more details on its definition).

Main results

Evaluating distortion correction is a difficult task since the non distorted image does not exist. At best, one can compare distortion correction results with a known non distorted image of another modality, although in that case direct comparison of the intensities is not possible. We have therefore evaluated our algorithm on phantom and in vivo data. Results on the phantom are not displayed here but available in [Hédouin et al. 2017] and demonstrate state-of-the-art results of our algorithm with some results better than TOPUP [Andersson et al. 2003] on some regions. In vivo results are presented in the following sections. They relied on a set of five images from control subjects acquired on a Siemens 3T scanner (images size: 128×128×60, resolution: 2×2×2 mm 3 ) with a total of 30 directions at b-value 1000 s.mm -2 . The images were acquired with 4 different PED: anterior-posterior, posterior-anterior, left-right and right-left. In both evaluations, the b 0 images were used for computing the correction which was then further applied to all DW images.

Visual evaluation

We have first evaluated the correction provided by our algorithm through the comparison of the corrected images with conventional anatomical images (3D T 1weighted images acquired in the same session as the control subjects). These results are illustrated in Figure 2 We then compared quantitatively our approach to two state-of-theart methods: [Voss et al. 2006] image lines registration approach, and TOPUP [Andersson et al. 2003] (available as part of the FSL package). We performed the following experiments on an Intel Xeon 2.5 Ghz computer on 20 cores. The mean computation time per subject is very short (about 5 s) for Voss et al. algorithm, 170 s for our algorithm and 500 s for TOPUP. Unlike TOPUP, BM is multi-threaded, allowing a faster computation time useful in the clinic.

To obtain a quantitative evaluation of the quality of the corrected images, we computed (for each subject) the distortion correction from the two pairs of EPIs (left-right/right-left, anterior-posterior/posterior-anterior). The idea behind this evaluation was that if the correction performs well, then the two corrected images should be the same (up to some additive noise). We have thus computed after correction the average of the local correlations between the two corrected images. These results are reported for the two methods in Table 2

Modeling water diffusion

After artifacts correction and image improvement, the next step is to infer the 3D probability of water diffusion from the diffusion sensitized measurements, which describes the white matter microstructure of the brain. It has been widely established (see Section 1.4 of [Stamm 2013] for more details and references) that diffusion measurements are directly linked to the probability of water diffusion through the image formation model illustrated in Figure 2.4. In this image formation model, we consider that measurements are performed in q-space, i.e. along vectors q which are functions of the gradient application time

Image = S 0 • |FT [p x ] (q)| Magnetic field gradient Image = • |FT [ ] (q)| Diffusion profile
We measure this

Baseline image

We want to infer that δ, time between successive gradients ∆, gyromagnetic ratio of hydrogen γ H , unit gradient direction g and gradient strength G. These q-vectors are often separated into two quantities: the b-value b = γ 2 H G 2 δ 2 (∆δ/3) and the gradient direction g. This simplification is made since almost all models do not need the separate parts of b. The attenuation observed from the baseline image then comes from the Fourier transform of the diffusion probability in space. Estimating the water diffusion in space after a certain amount of time thus amounts to solving the inverse problem from this image formation model. The optimization being performed on the signal, the models are often described by their Fourier transforms (i.e. characteristic functions), which we do in the following, but the parameters also describe the PDF. Before performing the estimation, the diffusion model, also called diffusion PDF or Ensemble Average Propagator (EAP), needs to be defined. Depending on the application and quality of the EPI acquisition, several models of different complexities have been defined. While this is not the main focus of this chapter, let us recall the main (non-exhaustive) categories of models in the literature, before explaining and focusing more on DCM and our contributions to their estimation.

Diffusion tensor

Historically, the Diffusion Tensor (DT) [Basser et al. 1994] was the first model introduced going beyond a model parameterized by a single scalar (Apparent Diffusion Coefficient (ADC)). It assumes the water diffusion PDF p x follows a centered (zeromean) multivariate normal distribution characterized by its covariance matrix D (symmetric positive definite):

D =   d xx d xy d xz d xy d yy d yz d xz d yz d zz   (2.5)
This model is one of the simplest to represent anisotropic diffusion, yet meaningful as it is considering anisotropic diffusion inside each voxel (see Figure 2.5). For these reasons, it is also one of the most used in clinic. Such a description also provides straightforward parameters of the tissue microstructure through the eigen analysis of the tensor D. Usual metrics [Basser & Pierpaoli 1996] combine the eigenvalues of D to cover a range of diffusion properties such as ADC, Axial Diffusivity (AD), Radial Diffusivity (RD) explaining the amount of diffusion in the voxel in all or specific directions, and Fractional Anisotropy (FA) quantifying the anisotropy of the diffusion. Figure 2.5: Illustration of some diffusion tensors with different diffusivities and anisotropies. Images are courtesy of [Kindlmann 2004].

Due to their simplicity, a large number of studies [Filippi et al. 2001] have therefore quantified the changes in those metrics due to several factors: diseases, aging, etc. However, the DT model suffers from two major drawbacks [Mori 2007]:

• A single multivariate normal distribution by nature assumes diffusion is happening in a single, principal direction, in a plane of directions or in all directions (isotropic). However, such a representation cannot handle properly a voxel in which fiber bundles with different directions cross. In this case, the estimated tensor does not represent the true nature of the underlying microstructure, which in turn affects microstructure studies or tractography algorithms.

• Microstructure parameters extracted from the DT are entangled: scalar measures, in particular ADC and FA, group several properties of the white matter into a single scalar value, and are thus difficult to interpret. For example, an FA or ADC change may be caused by an edema (i.e. inflammation that as a consequence brings free water locally inside the voxel) or by a specific destruction of axons in a fiber bundle. Using only the DT does not allow to separate easily those two sources of diffusion changes.

Chapter 2. White matter microstructure from diffusion imaging

Orthogonal bases

Although not covered in this manuscript, the first options to go beyond the tensor have consisted in using orthogonal bases representations of the signal. Since they were historically primarily presented for tractography and for an improvement of fiber bundles directions estimation, they are often called Orientation Distribution Function (ODF) or EAP. In this section, we consider all those models which start from modified Spherical Harmonics (SH) bases and go beyond but always share a common property/assumption: the signal is represented as a weighted linear sum of functions forming an orthogonal basis i.e.

S(b

i , g i ) = S 0 N j=1 c j Y j (b i , g i ) (2.6)
where (b i , g i ) represents the couple of the b-value and gradient direction of the i-th diffusion volume, S 0 represents the non diffusion weighted signal value, Y = {Y 1 , . . . , Y j , . . . , Y N } forms an orthogonal (or possibly orthonormal) basis, and c j (j = {1 . . . N }) are the linear coefficients of each basis function.

Illustration of the first elements of the modified SH basis are illustrated in Fig. 2.6. One of the earliest examples of such a model is the one proposed by [Descoteaux et al. 2007] where Y is chosen to be a modified, real basis of SH functions on the sphere [Atkinson & Han 2012]. This model allowed the estimation of several crossing fiber bundles directions per voxel, although being limited to a single shell acquisition due to the initial choice of basis functions on the sphere. These bases were then used to develop new and more specific tractography algorithms [Descoteaux et al. 2009]. On the other hand, this basis only considers fiber orientations and therefore loses microstructure information.

Other works in this category moved beyond the spherical harmonics and consider orthogonal bases that are compatible with multiple shells acquisitions. Among them, [Descoteaux et al. 2011] proposed the Diffusion Propagator Imaging (DPI) that extend the SH basis to multiple shells with the Laplace equation. Other approaches include the Simple Harmonic Oscillator Based Reconstruction and Estimation (SHORE) basis proposed in [Özarslan et al. 2013a] or Spherical Polar Fourier (SPF) expansion proposed by [Assemlal et al. 2009]. Finally, one may also consider R 3 as the basis and the fact that signal formation is directly related to the Fourier transform of the PDF of water displacement (see Fig. 2.4): this strategy, called Diffusion Spectrum Imaging (DSI) [Wedeen et al. 2005], thus takes the inverse Fourier transform of a large number of signals in q-space to get the PDF. While being model free, this last option however requires a very large number of signals to be robust which makes it difficult to use in clinical practice yet.

All these approaches have shown great interest for modeling diffusion, especially since the orthogonality of their bases makes their estimation rather simple and efficient (linear least squares problem to solve). Some rotationally invariant scalar properties of the models have been devised, e.g. return to origin probability [Özarslan et al. 2013b], however they do not convey a direct microstructure related property of the underlying tissues. Some work therefore remains to be done for making these models easily interpretable for clinical use.

Diffusion compartment models

Another class of approaches to recover the heterogeneous microstructure inside a voxel relaxes the constraint of a model being represented using an orthogonal basis. Rather, models in this class assume a given voxel is composed of different tissues (glial cells, water...) or fiber bundles with different orientations. Each of these groups of tissues, hereafter called compartments, are then assumed to have different properties, called microstructure properties, that characterize their water diffusion and thus their signal decay when a specific gradient is applied. Each voxel signal is thus assumed to be represented as a DCM, i.e. a weighted sum of compartments:

S(b i , g i ) = S 0 N j=1 w j ϕ j (b i , g i ) (2.7)
where ϕ j represents the signal decay of the j-th compartment (connected to a specific tissue type or fiber bundle with a specific direction), and w j represents the weight of the j-th compartment. An example of DCM is illustrated in Fig. 2.7. Again, while very similar to Eq. (2.6), this signal formation equation includes a fundamental difference in that every ϕ j explicitly represents a tissue type or main orientation in the voxel. This, with adapted models for each compartment, allows for the direct characterization of microstructure properties and thus potential changes due to development or pathologies, for each tissue at the sub-voxel level. For this reason, using DCM sounds promising for disease study, able to solve both problems of the simple DT model i.e. disentangling of parameters of diffusion change and crossing fibers resolution. Due to these advantages for disease characterization, developments have been made in this category, making it an active field of research. As for orthogonal bases, the first problem to tackle for DCM is the definition of the models assumed in each compartment i.e. for each tissue type or fiber bundle. This is a field of research in itself and it would go beyond the scope of this document to review all of them. Recent articles give an insight into quite a lot of the models that may be assumed for each compartment [START_REF] Panagiotaki | Compartment models of the diffusion MR signal in brain white matter: a taxonomy and comparison[END_REF]. For the sake of a better comprehension of future sections and chapters though, here are some typical compartments that may be encountered.

The first class of compartments is composed of the so-called isotropic compartments in that they represent water diffusing equally in all directions. The simplest and best example is the isotropic 3D Gaussian distribution parameterized by a single scalar diffusivity d. These compartments are usually used to represent 1-water diffusing freely around cells while being sufficiently far away from them not to be affected by their presence ; or 2-water diffusing inside cells whose shape can be globally considered spherical (e.g. neuron cellular part, glial cells, etc.). The second class of compartments is composed of the anisotropic compartments i.e. compartments associated to water diffusing principally along a given main orientation. These are typically compartments aligned with a given fiber bundle going through the voxel, up to a number of three anisotropic compartments: even though any number may be used, it is indeed commonly accepted that at most three fiber bundles cross in a brain voxel.

A very large number of anisotropic compartments have been proposed in the literature. Among the most known (non exhaustive list), multi-tensor models [Scherrer & Warfield 2012] use a classical tensor for each anisotropic compartment. Other models assume two subparts for a single anisotropic compartment to better explain the non Gaussianity of the signal decay for large b-values. These two subparts often correspond to the intra-axonal and extra-axonal parts of the diffusion for the fiber bundle. Two examples of these anisotropic compartments are the Neurite Orientation Dispersion and Density Imaging (NODDI) model [START_REF] Zhang | [END_REF] or the DDI model [START_REF] Stamm | [END_REF].

The previous advantages however come at the cost of a much more difficult non linear estimation problem, mainly due to the non orthogonality of the functions ϕ j , that remains an open problem. We have recently deeply studied this topic and present it in the next sections.

Estimation of diffusion compartment models

Estimation problem formulation

Model estimation follows usually the same framework for any diffusion model, aiming at minimizing the discrepancy between the model generated signals S(b i , g i ) for each b-value / gradient pair and the corresponding measured signals S i , which are assumed to follow a specific noise model. While some approaches have studied diffusion model estimation under the assumption of Rician noise [Fillard et al. 2007] or other noise models [Stamm et al. 2014b], most estimation frameworks assume Gaussian noise. In other words, we assume that for a given pair (b i , g i ) and S 0 signal, the signal measured is equal to the following:

S i = S(b i , g i ) + ε (2.8)
where ε follows a zero-mean Gaussian with variance σ 2 . In this specific case, the problem of estimation is formulated as a least squares system to be minimized over the parameters θ of the model and S 0 :

{ θ, Ŝ0 } = arg min θ,S 0 M i=1 (S i -S(b i , g i )) 2 (2.9)
where M is the number of DW images acquired, θ is the set of model parameters, θ and Ŝ0 are the optimal values of the parameters. Such a problem is typically minimized using different approaches depending on the model and the possibility to obtain derivatives of the cost function. For example, [START_REF] Panagiotaki | Compartment models of the diffusion MR signal in brain white matter: a taxonomy and comparison[END_REF] utilized a Levenberg-Marquardt optimization algorithm [Levenberg 1944] for this least squares problem for different DCM. Other approaches [Scherrer & Warfield 2012] used the BOBYQA optimization strategy [Powell 2009] to avoid computing too complex derivatives of the cost function. All these methods share the problem of trying to optimize a non linear least squares problem and are thus very sensitive to various artifacts in the acquisitions. In addition, DCM estimation is a very slow estimation problem due to the number of parameters to be estimated and the complexity of the cost function.

We have proposed a new approach [START_REF] Stamm | [END_REF], Commowick et al. 2016] towards 1-the simplification of the cost function and therefore a faster implementation of DCM estimation, and 2-a complete maximum-likelihood estimation framework that has, among other good properties, the ability of computing indirect values of the confidence in the estimated values, through the computation of the noise variance σ 2 (the larger σ 2 , the lower the estimation quality). Instead of the regular least squares formulation in Eq. ( 2.9), we formulate the estimation as the following maximum-likelihood estimation under a Gaussian noise assumption (see Appendix A):

{ θ, Ŝ0 , τ 2 } = arg min θ,S 0 ,τ 2 M 2 log τ 2 2π - τ 2 2 M i=1 (S i -S(b i , g i )) 2 (2.10)
where τ 2 = 1/σ 2 is the inverse of the Gaussian white noise variance on the input signals. Following the reasoning in Appendix A, this problem actually resorts back to solving the problem in Eq. ( 2.9) but additionally provides us with an estimate of the "local noise" as τ 2 which in fact contains both noise variance and errors coming from model unsuitability to the observed data.

Variable projection solution

The maximum-likelihood problem in Eq. ( 2.10) has the particularity of having some of its variables linear in the system (namely the weights of the individual compartments of the DCM). We have therefore developed in [START_REF] Stamm | [END_REF]] a variable projection solution [Golub & Pereyra 1973] to the optimization problem enabling fast optimization of the Maximum Likelihood (ML) estimator of DCM. In more details, our framework for DCM estimation considers a set of parameters composed of three subsets:

• parameters independent of the model: base signal S 0 , noise variance parameter of the ML estimation σ 2

• weights of the individual compartments of the DCM: w = {w 1 , . . . , w N }

• parameters θ j of the individual compartments of the DCM forming the set of parameters θ

Thanks to variable projection, we have developed a framework that, from the estimation of the non linear parameters θ alone, provides a complete estimation of all the aforementioned parameters. Moreover, we developed in Appendix A a simple derivative scheme over both the individual residuals and the cost function itself in Eq. ( 2.10), that provides: 1-a Levenberg-Marquardt optimization to provide robust estimation, and 2-a generic gradient-based estimation for any compartment type in the DCM where only the Jacobian of the individual compartments (component-wise derivatives D F of matrix F in Eq. (A.12)) with respect to their parameters need to be defined for the framework to be adapted to a new model.

Example: practical implementation for multi-tensor estimation For the sake of clarity and further experiments let us, without loss of generality on other models, focus on a specific model for estimation: the multi-tensor model. We consider here a DCM made of multiple tensors, each represented by the following individual ϕ j :

ϕ j (b i , g i |θ j ) = exp(-b i g T i D j g i ) (2.11)
where the parameters θ j of the individual compartment are all expressed in the tensor D j . A tensor has a total of six parameters that can be expressed in very different ways. For practical reasons of setting reasonable bounds to the parameters and get tractable derivatives, we have chosen to parameterize each D j by a main unit direction of diffusion e j,1 (two parameters in spherical coordinates), a supplementary Euler angle a j ∈ [0, 2π] in the orthogonal plane to e j,1 giving rise to the second eigenvector e j,2 , and three positive additive eigenvalues parts d j,1 , d j,2 , d j,3 such that eigenvalues of the tensor are ordered:

λ j,1 = d j,1 + λ j,2 , λ j,2 = d j,2 + λ j,3 , λ j,3 = d j,3 .
With these parameters, D j is parameterized as:

D j = d j,3 I 3 + d j,2 e j,2 e T j,2 + d j,1 e j,1 e T j,1 (2.12)
where I 3 denotes the 3D identity matrix. In our variable projection setting, the ϕ j (b i , g i |θ j ) constitute the F i,j elements of F. Getting the derivatives of the ML estimation formulation after variable projection thus only requires the derivation of D j against its parameters which can be readily obtained by differentiating Eq. (2.12). Additionally, since parameters are bounded, we applied Panagiotaki et al. strategy [START_REF] Panagiotaki | Compartment models of the diffusion MR signal in brain white matter: a taxonomy and comparison[END_REF] to unbound them for estimation with the Levenberg-Marquardt algorithm.

Main results

We have evaluated this estimation framework against synthetic and control subjects datasets. I provide here a short summary of the main results, however more details about these experiments and the datasets used can be found in [START_REF] Stamm | [END_REF]].

Evaluation on synthetic datasets

We have first performed an evaluation on a synthetic dataset (see Fig. 2.8) where we simulated multi-tensor models with:

• two isotropic compartments (one for free water, one for isotropically restricted water i.e. water inside spherical cells). No parameters needed to be estimated for these compartments apart from their weights.

• different numbers of anisotropic tensors at each voxel (from zero to three compartments) with different orientation configurations

• varying proportions for each compartment from one voxel to the other With these synthetic data defined, we then simulated DWI signals using the HCP gradient scheme [Essen et al. 2013] and added Rician noise (25 dB), after which we evaluated the estimation computation time with respect to the convergence to the right solution. We compared for this task different optimization algorithms, and different frameworks for optimization:

• Method A1: our framework: variable projection, Levenberg-Marquardt optimization [Levenberg 1944] with analytical derivatives

• Method A2: variable projection, Levenberg-Marquardt optimization with numerical derivatives

• Method A3: variable projection, globally convergent conservative convex separable approximation (CCSA) optimization [Svanberg 2002] with analytical derivatives

• Method A4: variable projection, derivative-free bounded optimization by quadratic approximations (BOBYQA) algorithm [Powell 2009]

• Method B: joint weights and DCM parameters estimation as proposed in [Scherrer & Warfield 2012] using the BOBYQA algorithm

We evaluated the correct evaluation of the models against the ground truth for each of these methods with two metrics: 1-mean square error of weights, 2mean square error of all D j . Each metric was computed separately for a given number of anisotropic compartments in the voxels (from 0 to 3). The evolution of these metrics for decreasing relative stopping criteria (from 1.0 × 10 -1 to 1.0 × 10 -13 ) is plotted against their computation time in Fig. 2.9. The main results confirm first a much accelerated estimation with the variable projection framework with the analytical derivatives and Levenberg-Marquardt optimization (method A1), while other gradient-based methods show longer computation times and method B shows computation times increasing much faster with decreasing tolerance levels. For zero to two fibers areas, method A1 is also outperforming all other methods, especially state-of-the-art method B. Method A1 indeed reaches its lowest level for higher tolerances and reaches lower mean square errors than other methods. All methods however seem more sensitive to initialization and convergence tolerances when considering three compartments which is likely due to much more complex estimation when considering three fibers.

8

* * * We have evaluated visually the ability of methods A1 and B to obtain good estimation results of multi-tensor models in a reasonable computation time. We thus optimized the parameters of estimation for both methods to get them to run in 30 seconds over a crop inside the corpus callosum and another one in the centrum semi-ovale (known for its fiber crossings). We report the results in Fig. 2.10. Visual inspection of these results further shows the ability of method A1 to provide more spatially coherent estimates of the multi-tensor models with less artifacts, mostly visible in three fascicle areas (inside the centrum semi-ovale).

Conclusion and perspectives

Artifacts correction in diffusion weighted imaging

We have conducted research to enable fast and accurate distortion correction of diffusion images. This work is crucial as it is a necessary first step for further accurate processing and studies of dMRI in a clinical context. Thanks to our work with Renaud Hédouin, we have presented a method that generalizes block-matching registration and provides state-of-the-art distortion correction or even better correction on phantom experiments (see [START_REF] Hédouin | [END_REF]). We have defined priors on transformations between local regions by constraining their formulation to match the distortion model in EPI. This is crucial as a priori constraints lead to more robust algorithms, contrarily to a posteriori regularization. Moreover, this is the first use of a symmetric setting in BM registration. This framework is applicable to other traditional registration problems, including atlas construction where a groupwise atlas construction can be derived from this symmetric registration.

Although less perspective can be directly seen from this work compared to other chapters, many options remain to extend and further evaluate this research. First of all, we have considered the most general case for our distortion correction (i.e. no anatomical non distorted image). In fact, we could definitely adapt our method to handle this kind of intermediate image. This would require extending the symmetric registration framework and change the similarity metrics used in the algorithm. We performed preliminary studies with Renaud Hédouin on that topic (see his PhD [START_REF] Hédouin | [END_REF]], Section 4.6), showing that it is doable but work remains to be done on the optimal similarity metric to define between three images. In addition to this direct extension, we could also extend this framework to treat different sources of distortion, namely Eddy currents [Mangin et al. 2002], by using similar local transformation definitions, this time to register individual gradient volumes.

Finally, evaluation of distortion correction has proven to be a difficult topic. For this reason, most articles tend to evaluate only visually their results which is not satisfactory. There is therefore a need for proposing new evaluation schemes, most probably based on phantom experiments or simulations, to properly evaluate differences in algorithms performance. Designing and using those simulations and phantoms could as well be future works, that could lead to a future challenge on distortion correction which would be very valuable to the field.

Diffusion models estimation

Model estimation for diffusion modeling is a critical step for further processing and statistical analysis of diseases. We have chosen over the past years to focus on diffusion compartment models as they provide a great, intuitive way of modeling the white matter microstructure. To overcome the challenges of DCM estimation, we have proposed a new ML estimation framework that allows for faster model estimation thanks to a variable projection formulation. Interestingly, this framework is versatile as it may be used virtually for any combination of compartment models, either using Levenberg-Marquardt estimation if the derivatives of the model can be computed or using gradient-free optimization otherwise.

While I decided to focus on this aspect of estimation in this manuscript, other problems arise very fast when considering the estimation of DCM. One of them concerns the determination of the optimal number of compartments in the model. We have proposed a new algorithm as a potential solution to this task [Stamm et al. 2014a], relying on the Akaike information criterion to determine the model over-fitting of the data. Validation is also an open problem. We have used here only a relatively simple simulation. No direct evaluation may be done on real life data. One way to further evaluate the model would be to adapt recent phantoms [Daducci et al. 2014, Caruyer et al. 2014] to simulate realistic multicompartment model data on several shells and evaluate our methods on them.

While the algorithms we proposed are a great step towards the use of DCM in a clinical context, some problems remain especially when a small number of gradient directions or b-values is available. Another critical point for model interpretation is the potential of some compartments to take the place of some others. For example, in a model with an isotropic free water compartment and several tensor compartments, nothing prevents one of the tensor compartments to take part or all of the isotropic weight by changing accordingly its parameters. This indetermination, coming from the fact that DCMs are not orthogonal bases, could lead to difficult interpretation of the resulting models as the free water weights thus cannot be compared directly between patients, especially when a large number of compartments are present in the voxel. Including priors on the compartments so that they cannot take over each other would therefore be an important step for the interpretability of the DCMs estimated. Ensuring a good estimation of the models when only clinical, low angular resolution data is available remains also a critical topic of research. We have in the past presented compartment models towards this objective [START_REF] Stamm | [END_REF]], but the careful optimization of estimation constraints remains an open problem to get a robust estimation as well as enough degrees of freedom to precisely evaluate the white matter microstructure.

Estimation of DCMs remains a long process which will need to be improved to be compatible with real time clinical use (one of the strengths of the DT model).

One solution to explore for this is dictionary-based initialization, inspired e.g. from [Yap et al. 2016]. Finally, DCM estimation is only a step towards our goal of getting more specific measures of white matter microstructure degeneracy. Tools are now needed to perform the processing (registration, interpolation) and statistical analysis of these models. I will cover some of these topics in Chapter 4.

Chapter 3

Multiple compartments T 2 relaxometry: insights into myelin

Quantitative relaxation times from MRI

Relaxometry gathers a set of acquisition methods whose aim is to measure, from a set of MRI signals, the relaxation times of tissues. Those relaxation times are at the basis of all contrast images (T 1 -weighted, T 2 -weighted...) and are quantitative i.e. these measurements theoretically do not change (at least much less than weighted images contrasts) depending on the scanner (apart from a change of B 0 field strength). This, plus the fact that relaxation times have a clear physical meaning, makes their measurement very interesting in practice. Two complementary communities have studied this topic: the acquisition community and the signal processing community. A large review of these aspects is proposed by [Tofts 2004].

Three main relaxation times are commonly measured through MRI (illustrated in Fig. 3.1): longitudinal relaxation time (T 1 ), and transverse relaxation times (T 2 and T *

2 ). Put shortly, T 1 corresponds to the time, after the resonance has stopped, at which the magnetization along the longitudinal axis (i.e. the nominal magnetic field B 0 direction) has reverted back to 63% of its original value before resonance. T 2 on its side concerns magnetization in the transverse plane. When resonance happens, all spins are put in phase in the transverse plane. When the excitation stops, these spins are going to dephase and therefore decrease transverse magnetization. T 2 corresponds to the time at which the magnetization would reach 37% of its original value after excitation stopped. T 2 corresponds to the ideal case if the acquisition was only influenced by the tissues imaged. For several reasons including inhomogeneities in the magnetic field, dephasing of the spins happens faster than the true T 2 value. Therefore, T * 2 is defined as the actual observed time at which the magnetization reaches 37% of its original value. In both cases, dephasing happens much faster than longitudinal relaxation and we therefore have the relationship that

T 1 ≥ T 2 ≥ T * 2 .
(a) (b) All of these constants have great interest for clinical purposes as tissues will have different T 1 , T 2 and T * 2 values depending on their composition. Therefore tissue changes due to pathologies are reflected in these three numbers. T * 2 is for example sensitive to changes in magnetic susceptibility (which results in changes of the B 0 field inhomogeneity). Its study gave rise to quantitative modalities such as Quantitative Susceptibility Mapping (QSM) [START_REF] Young | [END_REF], Wang & Liu 2014] which, although not the topic here, is used for many diseases e.g. MS where changes in the basal ganglia was observed early [Langkammer et al. 2013] or Alzheimer's disease where iron deposition is a known effect of the disease [Acosta- Cabronero et al. 2013].

The longitudinal relaxation time T 1 has also been widely studied in the literature. Techniques have been developed to obtain quantitative T 1 images first from gold standard techniques such as inversion recovery or saturation recovery [Crawley & Henkelman 1988]. Then dedicated acquisitions were proposed that are compatible with clinical acquisitions. Among many others, popular methods to measure T 1 include DESPOT1 [Deoni et al. 2005] or the MP2RAGE sequence [Marques et al. 2010], both relying on the acquisition of T 1 -weighted images with different flip angles. This type of clinically compatible acquisitions lead to a large range of applications for quantitative T 1 measurements. For example in MS, patients were found to have significantly higher T 1 relaxation values throughout the white matter [Vrenken et al. 2006b, Vrenken et al. 2006a] showing a disease activity, even in these apparently normal parts of the brain on conventional imaging.

The T 2 relaxation time has perhaps seen most of the developments in the past years. It provides, usually after non linear estimation from a sequence of multi-echo spin echo images (e.g. Carr-Purcell-Meiboom-Gill (CPMG) sequence [START_REF] Carr | [END_REF], Meiboom & Gill 1958]), complementary information to T 1 measurements highlighting changes in tissue microstructure. Some examples of the direct application of T 2 measurement include hippocampus study and the relationship between T 2 measurements and abnormalities leading to epilepsy [Sumar et al. 2011, Pell et al. 2004] or the study of MS lesions, normal appearing white matter and their evolution over time in MS [Combès et al. 2016, Kerbrat et al. 2017]. Other applications also include the study of pediatric T 2 relaxation times evolution depending on the brain region [Leppert et al. 2009].

In addition to single T 2 relaxation time estimation, more and more teams have considered the fact that each voxel in the brain is composed of several tissues with different T 2 values. As illustrated in [MacKay & Laule 2016] and shown in Figure 3.2, a voxel, due to its relatively large volume compared to the average cell size, is composed of a set of tissues each with different T 2 values: myelin has a short T 2 value in between 10 and 40 ms, intra-and extra-cellular matter (gray matter cells, axons...) have medium T 2 values around 100 ms, and free water far from any cellular structure has a large T 2 value around 2000 ms. All of these tissues and their proportions (see Fig. 3.2.c) are interesting, however one is crucial when studying the status or evolution of many brain diseases: the myelin. Myelin is indeed responsible for the fast transmission of the signal along axons and thus all brain function and other body functions depend on its integrity. A vast part of the literature of T 2 relaxometry has therefore looked at the possibility to estimate Myelin Water Fraction (MWF) from T 2 relaxometry sequences.

A large body of literature has studied the MWF in different diseases (MS [Laule et al. 2004], chronic stroke [START_REF] Borich | [END_REF], autism [Deoni et al. 2015]...).

Since it may very well help in solving the clinical-radiological paradox for MS, we have also studied this modality in depth. Two major issues however need to be tackled to enable robust quantitative markers of white matter microstructure in relaxometry:

• How to model the T 2 distribution in a voxel? Two different approaches fight each other on this topic: model-based where a priori constraints are set and model free where a posteriori regularization ensures estimation robustness. We debate on the current studies with both approaches in Section 3.2 and then go on with the second major problem.

• How to robustly estimate an Relaxometry Compartment Model (RCM)? Data acquired in relaxometry sequences, contrarily to diffusion in Chapter 2, are indirect measurements of the fractions of each tissue. This renders the estimation complex. We debate this aspect and present two approaches for RCM estimation in Sections 3.3 and 3.4.

Modeling multiple T 2 compartments

Many algorithms have been developed to estimate, from relaxometry sequences, either the MWF or the complete tissue microstructure at each voxel. Among them, methods using multi-echo relaxometry sequences are the most common. These multi-echo sequences consist in the acquisition of a series of 3D volumes, each for a different echo time e i and with a constant repetition time T R . Going further into specific acquisitions, the CPMG sequence acquires many images separated by a fixed echo spacing i.e. e i = i∆ TE . In its most general general form, the signal obtained at a given voxel at the i-th echo (i ≥ 1) is computed as follows:

S i = S 0 N j=1 w j R + p j (t)A(t, e i )dt (3.1)
where S 0 denotes the baseline signal if no attenuation was present (proportional to proton density), j denotes a j-th tissue category: p j (t) is its PDF, w j its weight ( j w j = 1). A(t, e i ) is the attenuation function that will be applied for a given relaxation time t at the echo time e i . Let us first discuss about p j and the number of compartments N . The choice of these two components has lead to two large families of estimation problems, although they can be summed up by the same equation.

If we choose in Eq. ( 3.1) a degenerate PDF for p j (t) as being a Dirac function δ(t j ) (i.e. null everywhere except at t j ), the integral simplifies itself and we get:

S i = S 0 N j=1 w j A(t j , e i ) (3.2)
Combining this choice with a large number of components N spread over the whole T 2 spectrum, we obtain the so-called multi-component model for the relaxometry signal. This approach has been used in many works including [Whittall & MacKay 1989, Prasloski et al. 2012, Layton et al. 2013, Dingwall et al. 2016]. All these algorithms have in common that they fix a large number of Diracs along the T 2 spectrum and estimate the weight of each of the pikes, usually through a non-negative least squares method. MWF is then usually obtained by considering the sum of all peak weights whose T 2 value is below a threshold compatible with myelin e.g. 50 ms. However, while simple in appearance, the estimation of the weights for so many peaks is highly under-determined and some regularization is necessary. Different regularizations in the literature include Tikhonov [Whittall & MacKay 1989], non-local [Yoo & Tam 2013] or spatial regularization [START_REF] Hwang | [END_REF], Raj et al. 2014].

The previous multi-component model imposes a regularization a posteriori of the obtained peak weights. Instead, another solution is the Relaxometry Compartment Model (RCM) which assumes a small number of compartments N , typically the three aforementioned compartments (short T 2 , medium T 2 and high T 2 water). Each of these compartments is now assumed to have p j accounting properly for the distribution of the class a priori and thus removes the need for regularization in the estimation problem. Taking an example, one may take a Gaussian PDF for each p j [Melbourne et al. 2013, Chatterjee et al. 2017a] with specific parameters modeling the typical aspect of each class as described in the literature [Laule et al. 2007]. For this class of algorithms, the problem therefore moves to the modeling and optimization strategy: choice of the right PDF [START_REF] Akhondi-Asl | [END_REF], Chatterjee et al. 2017b, Akhondi-Asl et al. 2014], choice to optimize or not the PDFs parameters [Layton et al. 2013], and how to perform this more complex estimation [START_REF] Akhondi-Asl | [END_REF].

Because of its a priori modeling properties and after investigation of the multi-component techniques through the internship of Lucas Soustelle [Soustelle et al. 2015], we have gradually chosen over the past years to investigate multi-compartment techniques (hereafter named RCM) with the PhD of Sudhanya Chatterjee.

In addition to these modeling and estimation aspects, one term has been discarded so far and has also been the topic of several research works: the attenuation term A(t, e i ). In perfect acquisition conditions, this term is a pure exponential term according to Bloch equations [Tofts 2004]:

A(t, e i ) = exp(- e i t ) (3.3)
Instead, CPMG sequences are subject to imperfect refocusing due to B 1 inhomogeneities [START_REF] Crawley | [END_REF]. This leads to stimulated echoes in the acquired signals starting from the second image. To handle this problem, the Extended Phase Graph (EPG) algorithm was proposed to model the attenuation over multiple echoes with imperfect refocusing [Layton et al. 2013, Prasloski et al. 2012]:

A(t, e i ) = EPG(t, i, ∆ TE , T 1 , B 1 ) (3.4)
where B 1 is a scalar multiplicative factor that models the imperfect refocusing and has to be estimated from the data, T 1 is the longitudinal relaxation time at the current voxel (which can be estimated from a quantitative T 1 relaxometry sequence). This EPG attenuation and the pure exponential one are illustrated in Fig. 3.3 and demonstrate the ability of EPG to handle signal from stimulated echoes. For this reason, although a topic of research on its own, we will only consider in the following the EPG attenuation as it fits well the acquisition model.

Robust compartment models estimation

We have described a new method for the robust estimation of RCM [Chatterjee et al. 2017a], from clinical data. I present here a brief introduction to the method used and the main promising results we obtained on a preliminary study on MS patients.

Non-negative compartment weights estimation

We use in this work an RCM tuned for the use on clinical data, where time constraints lead to the acquisition of a restrained number of echoes for the T 2 relaxometry sequence. Typically, retrospective studies on diseases only consider seven to eleven echoes acquired with a CPMG sequence with an echo spacing around 10 ms. Additionally, even with a large number of echoes, [Layton et al. 2013] have highlighted the difficulty to estimate the parameters of T 2 distributions. We have therefore chosen in this work to consider a set of 3 compartments each with a PDF p j with fixed parameters tuned to model the three compartments of water that may encountered in a voxel: short, medimum and high T 2 . This approach is in fact independent of the PDF chosen for p j and we now detail the estimation procedure for any PDF choice. At a given voxel, the unknowns of the estimation problem in Eq. 3.1 are now the baseline signal S 0 , the B 1 inhomogeneity factor and the weights w = {w 1 , w 2 , w 3 }. We first choose to fuse S 0 and w into a vector of variables α, as in Appendix A, whose terms are only constrained to belong to R + . From this vector, S 0 is obtained as the sum of α, and w is obtained by dividing α by S 0 . We therefore get the following optimization problem:

α, B1 = arg min α,B 1 Ne i=1   y i - N j=1 α j R + p j (t)EPG(t, i, ∆ TE , T 1 , B 1 )dt   2 (3.5)
where N e is the number of echoes acquired. While solving this problem for both B 1 and α is complicated, solving for each variable independently is simple. We have thus chosen the following alternate optimization scheme:

• Fix B 1 and optimize over α: this is achieved by using a non-negative least squares algorithm [Lawson & Hanson 1995],

• With this updated α, optimize B 1 : this problem is non linear in nature and its derivatives are computationally expensive to compute. We therefore optimize it using the gradient-free BOBYQA algorithm [Powell 2009].

Main results

We have modeled each compartment distribution (p j ) using a Gaussian PDF, each centered around the typical values of T 2 : for the short T 2 (mean: 20 ms, standard deviation: 5 ms), medium T 2 (mean: 100 ms, standard deviation: 10 ms), and high T 2 compartments (mean: 2000 ms, standard deviation: 80 ms). We present in the following the main results we have obtained with this model first on numerical phantom simulations, then on a longitudinal study on MS patients [Chatterjee et al. 2018a].

Evaluation on numerical simulations

We have performed many numerical evaluations that may be seen in Sudhanya Chatterjee's article [Chatterjee et al. 2018a] and PhD [Chatterjee 2018]. I present here only a representative evaluation we have performed on a numerical phantom illustrating the potential of the method. To evaluate the robustness of our estimation framework in different noise, B 1 and compartments configurations, we have built a numerical phantom of T 2 relaxometry sequences. We built the numerical phantom to be as close as possible to the real signal formation model. For each pixel in Figure 3.4, different configurations of compartment weights and B 1 value were selected to explore the range of values relaxometry signals may reach in a real case.

For each pixel, different T 2 curves (each of 32 echoes each spaced by ∆ TE = 10ms generated using EPG with the specified B 1 factor) were randomly sampled along the distribution and averaged to obtain realistic signals of combinations from molecules constituting the tissues. From these reference signal images and the known ground truth T 2 distribution, we simulated different levels of Gaussian noise to obtain a Signal to Noise Ratio (SNR) between 50 and 1000 (typical relaxometry images have a SNR between 50 and 100) and run our estimation algorithm to obtain the compartment fractions (the compartment PDF parameters being the same between the simulation and estimation algorithm). Validation of the robustness of the estimation method was computed as a relative mean square error between the true parameters and the estimated ones. This relative error is presented as a function of SNR in Figure 3.5.

Examining these results lead us to the conclusion that the proposed algorithm finds reliably the weights of each T 2 compartment. It is indeed visible on Fig. 3.5 that the relative mean square error does not go above 0.01 even at very low SNRs, indicating a good reliability of the estimates over noise in all configurations. For the same phantom, we also generated T 2 relaxometry signals using only 7 echoes each spaced by ∆ TE = 13.7ms, a number of echoes more acceptable for clinical acquisitions. We then estimated water fractions in each of the T 2 compartments and measured the relative absolute difference between the 7 and 32 echoes sequences. These results are presented in Figure 3.6. This figure demonstrates the similarity between the weights estimated for each compartment, showing the robustness of the approach to a drastic reduction of the number of echoes and the increase in their time spacing. From our experiments [Chatterjee 2018], the only significant differences between 32 and 7 echoes sequences are seen for a low SNR of 50, which is the lower bound of the noise we have encountered in our images where SNR is more around 75 or 100.

Longitudinal study on MS patients

Based on the experiments on simulated data and the robustness shown to a change of protocol, we have moved on to applying the developed technique to the analysis of MS patients, on a specific database acquired at Rennes University Hospital. This database consisted of 10 patients followed from their clinically isolated syndrome over a period of three years. During these three years, scans were performed every three months during the first year, then every six months during the second year, and finally a scan after three years. At each of these scan sessions, images acquired included a Gd enhanced T 1 -weighted image to highlight active lesions, a FLAIR and T 2 -weighted image to delineate lesions, and a T 2 relaxometry sequence with 7 echoes spaced by ∆ TE = 13.7 ms. The resolution of the T 1 -weighted and FLAIR images were isotropic at 1×1×1 mm 3 , while the T 2 -weighted and relaxometry images were acquired on the axial plane with a resolution of 1.5×1.5×3 mm 3 . Additionally, an expert radiologist delineated on the first time point the lesions and the active parts of the lesions. This gave rise to two different sets: L-: MS lesion parts at the first time point that did not present an enhancement with Gd (i.e. inactive lesion regions), E+: MS lesion parts at the first time point that presented Gd enhancement (i.e. active lesion regions). This separation was made on the first time point and the delineated regions were then followed over the three years as shown in Figure 3 Based on this database, we have evaluated the ability of compartment fractions to distinguish these two kinds of lesion parts in terms of their evolutions. An example of the evolution of these values for a given lesion is presented in Figure 3.8. This figure illustrates the potential of water fractions to fully characterize the evolution of edema and myelin recovery in lesions. After an initial state where the lesion shows no visual sign of myelin and a clear sign of edema, the recovery becomes more and more pronounced over the scans, showing a gradual, partial recovery both in terms of high and short T 2 water fractions. We have continued the investigation over all patients and all regions of E+ and L-in the database. To do so, we have evaluated, for both groups (E+ and L-), the average short, medium and large T 2 water fractions over the period of 36 months. Parts of the results of this study are presented in Figure 3.9, a more detailed view may be found in [Chatterjee et al. 2018a]. Overall, these results show different trends of evolution for E+ and L-lesions, coherent with assumptions from the clinic. At the first time point, E+ lesions clearly have less myelin content (short T 2 water fraction) and more edema (large T 2 water fraction). Additional experiments on successive time point differences between lesion groups showed that E+ and L-lesions differ in their rate of recovery: L-lesions tend to recover faster their large T 2 water fraction than active lesions (E+), indicating more severity of the active lesions. At the end of the 36 months, both groups come back to similar trends of values. Medium T 2 water fractions do not change over time between E+ and L-. However, medium T 2 is more of an intermediate class gathering many different tissues and the contents of that compartment could vary over time without having its fraction change. Overall, these results are very encouraging on the use of relaxometry for the distinction of aggressiveness of the lesions in the patient, thereby providing robust and specific markers of the disease.

Towards compartment parameters estimation

While the previous method has demonstrated its robustness and potential in clinical studies to characterize lesions, it leaves behind a part of the problem: the estimation of the compartments parameters describing each T 2 distribution and therefore intrinsic tissue composition change in lesions. Doing that estimation is however very difficult as the optimization of these parameters is very sensitive to noise and other artifacts. Moreover, the distribution chosen in Section 3.3 for each compartment is a Gaussian, whose support is not on R + . On the positive side however, we can restrict ourselves to the estimation of only a subpart of the parameters and the estimation problem, as in Chapter 2, is linear in some parameters (compartment weights). We have therefore proposed in [Chatterjee et al. 2018c] a new approach towards PDF parameters estimation for RCM, very similar in spirit to DCM estimation.

Gamma relaxometry compartment model

Since a Gaussian PDF has a support on R, it is not the most adapted distribution to model tissue T 2 distributions (p j in Eq. (3.1)), especially for the short T 2 compartment whose mean T 2 lies relatively close to 0. We have therefore changed for a Gamma distribution that can fit well the expected shapes of the distributions and has a support on R + . Additionally, we have parameterized the Gamma distributions by their means and standard deviations (instead of the classical shape and scale parameters) for easier parameter setting and interpretation:

p j (t) = t -1+µ 2 j /σ 2 j σ 2 j µ j µ 2 j /σ 2 j Γ µ 2 j σ 2 j exp - µ j t σ 2 j (3.6)
with µ j : the mean of the Gamma PDF, and σ j : its standard deviation. Putting this back in Eq. (3.1) and using the same idea as in Section 3.3 for handling S 0 and weights, we then get an updated signal model whose parameters are α, B 1 ,

θ = {µ 1 , σ 1 , . . . , µ N , σ N }.

Maximum likelihood estimation framework

If we look at our signal formation model closer, an interesting thing is its similarity with the one in Chapter 2. While the individual compartment formulations are of course different, its structure is the same. Following the approach in Appendix A, we have then developed a maximum likelihood estimation framework. It amounts to solving the following least squares problem:

α, B1 , θ = arg min α,B 1 ,θ Ne i=1   y i - N j=1 α j R + p j (t)EPG(t, i, ∆ TE , T 1 , B 1 )dt   2 (3.7)
For the same reason as in Section 3.3 (complex derivatives with respect to B 1 ), we have again chosen an alternate optimization first fixing B 1 and estimating α and θ, then fixing α and θ and estimating B 1 . While the B 1 optimization remains similar as in the previous sections, we have taken advantage of the linearity of the α parameters by using a variable projection framework as presented in Appendix A simply by identifying the elements F i,j of the problem matrix F to the following:

F i,j = R + p j (t)EPG(t, i, ∆ TE , T 1 , B 1 )dt (3.8)
Only the derivatives with respect to nonlinear parameters need to be defined. We have performed this derivation [Chatterjee et al. 2018c] after a study of the cost function behavior as a function of its parameters and in the presence of noise. From this analysis we have demonstrated that, even though all derivations to solve the problem can be perfectly done, estimating all parameters for all p j was too unreliable and sensitive to noise. This probably comes from the fact that the observed signal is not a direct measurement but rather an integration over the compartment distribution. Our approach towards parameter estimation therefore considers the only parameter that was stable enough: µ 2 the mean of the medium T 2 compartment. After some calculations, we obtain the following expression for the derivative of the cost function against any of the mean parameters:

∂F i,j ∂µ k = R + ∂p j (t) ∂µ k EPG(t, i, ∆ TE , T 1 , B 1 )dt (3.9) ∂p j (t) ∂µ k = p j (t) µ j σ 2 j 2 log t µ j σ 2 j -2Ψ µ 2 j σ 2 j + 1 - t σ 2 j (3.10)
if k is equal to j, and 0 otherwise. In Eq. (3.10), Ψ is the digamma function.

Main results

We have evaluated this method on several aspects including two main ones: repeatability and evaluation of mean parameter variation in MS lesions.

Repeatability experiments

Repeatability experiments were performed to evaluate the robustness of parameter estimation, even when freeing the estimation of some of the individual PDF parameters. We have considered a test-retest acquisition of T 2 relaxometry data of 4 healthy controls. For each subject, the acquisition was performed twice moving the subject out of the scanner in between the two scans. The acquisition details were as follows: image size of 192 × 192, voxel resolution of 1.1×1.1×5 mm 3 , 32 echoes spaced by 9ms. On these datasets, fifteen white matter regions were marked (see Figure 3.10) on the first (test) scan and the values obtained for the different parameters were compared between the test and retest scans.

peatability of the proposed method. The age of the healthy controls was in the range of 26-32 years. 15 regions of interests (ROIs) were marked in the brain for each healthy control over which the test and retest values of the compartments' water fractions were compared. All the ROIs were marked for one case. The ROIs were then registered on the other cases using a rigid followed by an affine registration [Ourselin 2000[START_REF] Chatterjee | Multi-Compartment T2 Relaxometry Model Using Gamma Distribution Representations: A Framework for Quantitative Estimation of Brain Tissue Microstructures[END_REF]] to ensure that similar regions were analyzed for repeatability in all cases. An illustration of these ROI on a subject is shown in Fig. 6.5. Based on these regions, we have performed [Chatterjee et al. 2018a] a quantitative study of the repeatability (between the test and retest scans) of the compartment weight measurements through two techniques: Bland-Altman plots and linear regression on the measures. We illustrate in Figure 3.11 this second evaluation, whose results are very similar to the ones of the Bland-Altman plots. Those results highlight on each graph the deviation between the ideal identity regression and the observed one. On all parameters, we observe only a very small deviation between the two regressions which highlights the fact that the proposed algorithm is well repeatable on these datasets.

Multiple sclerosis patient case

In addition to this first evaluation, showing a stable estimation in terms of weights on test-retest data, we have then evaluated the added value of the proposed framework (especially of the estimated mean T 2 compartment parameter) on an MS patient.

We have studied a T 2 relaxometry MRI scan with multiple lesions with the following acquisition parameters: image size: 192×192 (single slice), voxel resolution: 1.1×1.1×5 mm 3 , 32 echoes each spaced by 9ms. We present in Figure 3.12 the obtained water fraction maps of the three T 2 compartments as well as the medium T 2 mean relaxation time as estimated for the patient. On this figure, several trends may be observed. First of all, as for the previous model with fixed parameters, a decrease of the short T 2 water fraction is observed indicative of demyelination in the lesions. This observation is coupled with a clear change of the medium T 2 compartment mean relaxation time: the values of this parameter are clearly larger inside the lesions when compared to normal appearing white matter. To further characterize lesions based on this last parameter, we have then explored for two lesions in Figure 3.12 the profiles of the PDF mean parameter as we cross the lesion in different directions. Three profiles are illustrated in Figure 3.13. On this figure, we can observe different profiles of variation depending on the direction of the profile or lesion. Profiles (a) and (c) indicate a lesion resulting from the fusion of two lesions (Fig. 3.13.c). For both of them, a clear increase of the mean relaxation time of medium T 2 is seen, indicative of a change in this compartment when going toward the center of the lesion. This is probably related to a recent lesion with more aggressiveness in the lesion center. On the contrary, profile (b) in Fig. 3.13 indicates a more homogeneous lesion with small or even no change at all when going towards the core of the lesion. These results are very interesting as they highlight some variations in the behavior of lesions, probably related to their aggressiveness (although this fact remains to be investigated).

Conclusion and perspectives

Myelin content inside brain tissues of a patient is a critical indicator of the status and evolution of many neurodegenerative diseases. Obtaining ways to quantify its presence or absence, at least indirectly, is therefore an important step towards a better understanding of these diseases. We have therefore proposed, particularly thanks to the PhD of Sudhanya Chatterjee but also to previous internships such as the one of Lucas Soustelle, new methodologies towards the robust estimation of multi-compartment models of T 2 relaxometry (RCM). These methods rely on two new frameworks for estimation: one very constrained but also very robust estimating only the three compartments' weights; the second one able to provide more information on the individual compartments parameters. Both approaches are or can be directly based on similar frameworks as for DCM using variable projection for an improved robustness.

We have demonstrated that both methods produce repeatable compartment weights on test-retest experiments and also on different acquisition sequences (variable echo spacing and number of echoes). These results indicate that both algorithms are very well usable and robust. We have so far applied these methods to a study of MS patients first at a single timepoint showing that both water fractions and mean PDF parameter change 1-with a gradient when going inside the lesion; and 2-when considering different lesions, in conformity with histological observations on different natures of lesions [START_REF] Lassmann | Heterogeneity of multiple sclerosis pathogenesis: implications for diagnosis and therapy[END_REF]. The second study we have performed concerns the longitudinal analysis of the evolution of lesions that also highlighted different patterns of evolution depending on the activity (in the sense of Gd presence) of the lesion at the starting timepoint: active lesions appear on this small sample to recover slower than inactive lesions, although those two groups are not different after some time.

One topic of future research will concern more in-depth studies of the true relation ship between tissues constituting the voxels and the obtained water fractions. This is especially true for water bound to myelin (the short T 2 compartment) that also includes other specific tissues in the brain or iron depositions. Methods to quantify myelin are also not very reproducible from one scanner to another especially if changing vendors. This does not impair the interest in myelin related measurements from T 2 relaxometry as they have proven very useful for studies of diseases. There is however a need for a more in depth study of the quantities measured by relaxometry and their relationship with true known tissues. Designing synthetic experiments, phantoms or coupling acquisitions with histology acquisitions for that would be very important. The latter option seems more and more doable as more repositories come out providing (for now only diffusion1 ) both MRI and histology acquisitions, that are very valuable for this task.

While the proposed methods provide robust results when considering only weights estimation, it is still very challenging to estimate the PDF parameters from the current acquisitions. Our studies with Sudhanya Chatterjee [Chatterjee 2018] have shown how estimating several parameters of the T 2 compartments is very hard and sensitive to noise and other artifacts. This problem of robust estimation of all parameters remains an open one that will need proper new methods to be resolved. Once resolved, our preliminary results on just one parameter suggest that we could imagine studies as for diffusion where the informative microstructure parameters are not only the compartment weights but a combination of weights and internal compartment changes (that cannot be seen nowadays).

On the path of our research, we have also seen how much diffusion imaging and relaxometry are complementary. First, diffusion imaging is blind to myelin due to a still too long echo time compared to the one of myelin. However, diffusion imaging provides information on the directionality of tissues and some of their microstructural properties that relaxometry is unable to provide. Combining the two modalities in joint evaluation frameworks will now be possible in a near future, and has in fact been preliminarily started as I will present in the following chapter. Those joint frameworks have also been started in the literature (e.g. gratio [Campbell et al. 2018]) although they do not fully exploit the measurements of both modalities but rather try to replicate histology measurements. We hope with these new evaluations to provide comprehensive analyses of directional, microstructural and myelin contents of the tissues. Additionally, this complementarity exists at the estimation level. Although not used at its best in the current versions of the RCM estimation methods, variable projection could be fully used for all methods presented in this chapter. One nice feature of this chapter and Chapter 2 is thus that they rely on exactly the same model structures (multiple compartments) and therefore same estimation frameworks. In a longer term, it would thus be interesting to study joint diffusion and relaxometry estimation to take advantage of all the information at once [Kim et al. 2017, Benjamini & Basser 2016, Canales Rodriguez et al. 2018]. This would allow for more robust estimation frameworks but would also require compatible sequences and updates to the estimation method.

Chapter 4

Quantitative image processing for disease study

Contents

Last but not least, this chapter explores our research on quantitative MRI processing, including interpolation, atlasing and combining modalities, to go towards a better comprehension of disease status and evolution. This work was conducted mainly with two PhD students: Sudhanya Chatterjee and Renaud Hédouin, in close collaboration with the Children's Hospital Boston (team of Simon Warfield). Among the papers that arose from this work, two are particularly detailed here: • Relaxometry and dMRI for active lesions detection: Sudhanya Chatterjee, Olivier Commowick, Onur Afacan, Simon Warfield and Christian Barillot.

Identification of Gadolinium contrast enhanced regions in MS lesions using brain tissue microstructure information obtained from diffusion and T2 relaxometry MRI. In MICCAI proceedings, pages 63 -71, 2018.

Interest of quantitative image processing tools

We have seen in previous chapters how diffusion imaging and relaxometry, among other quantitative modalities, are promising in characterizing diseases. These modalities indeed provide interpretable quantities that may inform the clinician on the disease status. When performing population studies, or when evaluating a patient against a set of healthy subjects, two options are generally available: the first one is to extract, from the diffusion model images, scalar maps that represent microstructure properties. These scalar maps can in turn be analyzed for example through their registration in a common reference frame, i.e. an atlas. Such analyses could greatly gain from the new models provided by quantitative images, allowing for better matches of the images and more interpretable parameters of the brain microstructure. In the end, quantitative analyses could lead to much improved specificity and interpretability for the clinician. Let us take, as an example, the case of diffusion imaging. There have been many studies on DT images processing that have shown great interest in using the full information of the tensor rather than scalar maps. First, many algorithms have demonstrated a much improved alignment of images when considering the full tensor compared to anatomical images [Ruiz-Alzola et al. 2002, Zhang et al. 2006, Yeo et al. 2009], the main reason being that the diffusion tensor depicts directional and microstructural information in the white matter where anatomical images provide uniform intensities. Following in this track, statistical comparison methods were proposed, generalizing several tests to tensor images [START_REF] Lepore | [END_REF], Whitcher et al. 2007] or performing longitudinal analyses [Grigis et al. 2012, Keihaninejad et al. 2013], again demonstrating better specificity in using the complete tensor information rather than part of it. All these developments were supported by great works on underlying processing frameworks on tensors, either defining mathematical operations in a Riemannian manifold [Pennec et al. 2006] or Lie group [Arsigny et al. 2006b] (allowing operations such as interpolation) or by studying their re-orientation after transformation (finite strain re-orientation [Ruiz-Alzola et al. 2002] or preservation of principal direction [Alexander et al. 2001]).

Based on this interest for diffusion model processing, several groups have been exploring processing methods on more complex models such as ODFs. Among those works, a Riemannian processing framework has been defined on ODFs and applied for the study of differences between populations [Goh et al. 2011]. Registration of ODF images was also explored by several works [Raffelt et al. 2011, Du et al. 2012]. However, DCM image processing remains yet to be really explored. Only a few papers study their processing and alignment [Taquet et al. 2014] and they remain often linked to a specific type of compartment model, preventing their use on the vast majority of models available in the literature [START_REF] Panagiotaki | Compartment models of the diffusion MR signal in brain white matter: a taxonomy and comparison[END_REF]. Such DCM image processing would however allow even more specific analyses and atlases, as suggested by previous studies on other models [Goodlett et al. 2006, Zhang et al. 2007]. I will therefore present some of our research on this topic in Section 4.2.

Towards DCM registration and atlasing 49

Alongside these developments, we have seen earlier that one of the interests of dMRI is to allow the analysis of microstructure changes along fiber tracts, especially along those that are known to be crucial in patient motor function or particular cognitive skills for example. Such studies on large groups could allow a better understanding of neurological diseases and thus help to resolve the aforementioned clinical-radiological paradox [Guttmann et al. 1995]. Some previous research, named tract-based spatial statistics (TBSS) already went in this direction [START_REF] Smith | [END_REF], Jbabdi et al. 2010], although using only the DT model which is too simple and whose parameters are too entangled, or very simple compartment models. Also they do not consider fiber extraction per se but rather voxel statistics on skeletons of the white matter architecture. With all these remarks in mind, we have decided to go towards a personalized evaluation of the patient disease burden at the compartment level along fiber tracts. I present these developments in Section 4.3.

In addition to the previous developments on dMRI, there is a great interest in the combination of information from different modalities to fully understand tissue characteristics in each area of the brain. A very good example of the complementarity of indicators from different quantitative images is the so-called g-ratio. Several papers [START_REF] Stikov | [END_REF], Campbell et al. 2018] have demonstrated the capability to compute, from the combination of diffusion imaging and myelin sensitive images, this number that corresponds in histology to the ratio of the axon plus myelin diameter and the axon diameter. This ratio, although subject to a debate of interpretations, suggests the interest of using relaxometry images and dMRI for getting the full picture of e.g. disease progression in a patient with myelin destruction. We have studied, and I present in Section 4.4, the joint processing and learning from these images for an also important project: the detection of Gadolinium (Gd) active lesions, without the need of Gd injection which has been shown recently to accumulate in brain tissues across time [Gulani et al. 2017].

Going towards diffusion compartment model images registration and atlasing

Registration of DCM images, like every registration algorithm, relies on several key parts such as the definition of the transformation being sought, how to optimize it, or the similarity metric defining how well the images are in correspondence. While these core components are underlying every registration algorithm, two additional key parts need to be defined both for the registration of oriented models (such as DCM images) and the constitution of atlases: 1-an interpolation / averaging scheme, and 2-a re-orientation strategy when applying a transformation. While the re-orientation strategy, extensively studied for tensors, can be directly applied compartment-wise to DCM, the interpolation / averaging scheme needs to be redefined. Such a scheme is of importance for example when applying a transformation to an image or when picking the diffusion model at a current position in a tractography algorithm, or even when computing an average DCM image in atlas construction. The second problem, not covered here, is the definition of similarity metrics between DCM images. I [resent in the following methods we have introduced to solve the first problem with a focus on being generic enough so that it can be applied to virtually any compartment model in the DCMs.

Diffusion compartment models interpolation and averaging

Combining images into an atlas or interpolating a position between several voxels requires the same operation. In the first case, a set of DCMs each weighted by 1/N , where N is the number of DCMs, has to be merged into a meaningful average DCM.

In the second case, interpolation provides the value at a given position by combining weighted DCMs from the neighboring voxel positions. In both cases, the average DCM has to be meaningful i.e. keeping as much information as possible from the input data, yet remaining simple enough to be computationally tractable. If we push further and use a simple idea, we can compute an average DCM from a set of weighted DCMs as the reweighted sum of all compartments that constitute the weighted input DCMs. This would keep the total information of the input data, but would also lead to non tractable models (e.g. averaging 100 DCMs for atlas construction, each with 4 compartments, would lead to 400 compartments). We have instead chosen with Renaud Hédouin to define the interpolation / averaging problem as a simplification problem, trying to keep as much information as possible while having a reduced number of compartments. Our idea, devised in [Hédouin et al. 2015], is similar to that of [Taquet et al. 2014, Taquet et al. 2015], yet being more generic in the compartments it can be applied too. Our simplification problem is shown in Our interpolation scheme is split into two parts. We first consider that at most one compartment of each isotropic compartment type is present in each of the input DCMs. For those, the simplification is therefore straightforward: only their parameters need to be averaged per compartment type and the weights computed, no further simplification is needed. The second part of our scheme concerns directional compartments. These compartments, assumed in our framework to all follow the same model (e.g. anisotropic tensor), are much more numerous (generally up to three for each of the input data). We therefore apply a simplification of these input compartments into a predefined number of compartments. This simplification is performed by clustering the compartments using fuzzy spectral clustering [START_REF] Ng | [END_REF] and then averaging each cluster into a single directional compartment. Weights are then recomputed from the input data to provide the final averaged DCM. This simplification approach is generic in nature. It only needs the definition of two parts for the application of the method to a given directional compartment type: 1-the definition of a distance between compartments to compute the initial affinity matrix for clustering, 2-the definition of compartment averaging into a summary compartment for cluster aggregation. We have defined in [Hédouin et al. 2015[START_REF] Hédouin | [END_REF]] these two key points and evaluated our scheme for the multitensor model where each directional compartment is an anisotropic tensor (using log-Euclidean distances and averaging schemes [Arsigny et al. 2006b]), and for the Diffusion Direction Imaging (DDI) model with several different metrics [START_REF] Stamm | [END_REF]].

Main results

We have evaluated our interpolation and averaging on several tasks and using different evaluation measures in [Hédouin et al. 2015[START_REF] Hédouin | [END_REF]]. I present a short summary of some of these results: 1-evaluation of two different metrics for interpolation of DCM images, 2-construction of an atlas from DCM images.

Evaluation of diffusion compartment model interpolation

We have first considered a simple experiment to evaluate our interpolation method. The diffusion images used in the experiment came from the HCP, where dMRI was acquired with a total of three b-value shells (from 1000 to 3000 s.mm -2 ) and 270 gradient directions, each volume being of size 145 × 174 × 145 and voxel size 1.25×1.25×1.25 mm 3 . From these images, we used the previously presented estimation algorithms (see Chapter 2) to estimate a DCM at each voxel with the following characteristics: one free water compartment, one isotropic restricted water compartment, and three directional compartments each modeled as a DDI compartment [START_REF] Stamm | [END_REF].

We have applied three successive rotations to the images, each having the same axis and an angle of 120 degrees. After these three rotations, the DCM image should come back to its original state and, if the interpolation is right, the difference between the two images should be as small as possible. Since the DDI compartment is complex, we have evaluated several options to compute distances and averages between them (four options in total). We illustrate the results of two of the best methods in Fig. This figure illustrates the differences arising from interpolation metrics. While both methods appear to provide good results compared to the original image (top line of Fig. 4.2), the rest of the figure highlights that in fact method 1 provides results closer to the reality in some aspects (and in fact smaller errors compared to the original image in terms of simulated signal from the models, see [Hédouin et al. 2015]). Going in more depth, orientations of the individual compartments are better with method 1, while their microstructure properties seem better with method 2. The choice of the distance used in clustering and of the averaging scheme of individual compartments is therefore crucial for obtaining correct interpolation results. The definition of what is correct is also very important: as [Taquet et al. 2015] also mentioned in their study, one question is whether we want to preserve microstructure properties or diffusion signal properties. This is an open problem that probably does not admit one solution, but rather the best option to take depends on the context.

Construction of a diffusion compartment model atlas

After evaluation of the interpolation scheme, we have then carried on with applying the technique to build an atlas of DCM images. To do so, we have followed [Guimond et al. 2000] atlas construction, adapting to use diffeomorphisms encoded as SVFs and the log-Euclidean framework [Arsigny et al. 2006a]. To compute the atlas, registration is needed. Since no DCM registration was available, we have chosen to build the atlas using the registration framework detailed in Appendix B, applied to DT images [Suarez et al. 2012] first (with an adapted similarity metric and re-orientation scheme [Alexander et al. 2001]). Then, as a post-processing step, we have applied the obtained transformations to the DCM images using our interpolation.

We have considered a database of 46 diffusion images with an image size of 128×128×55, voxel resolution of 2×2×2 mm 3 and 30 gradient directions with a b-value of 1000 s.mm -2 . From these, DCM images were estimated with one free water compartment and three DDI compartments at each voxel and used to build the atlas. The obtained atlas is illustrated in Figure 4.3. This atlas provides a clear distinction of crossing fibers and will be of great interest in future studies, for example in the ones presented in the next section.

Patient to population comparison of diffusion properties along white matter tracts

We have defined an interpolation method that permits the computation of DCM atlases as reference frames for further studies, for example comparing groups of patients and controls, or a single patient with respect to an atlas of control subjects.

We are now primarily interested in the latter option as it may provide, with the right tools and models, comprehensive and specific information on the status of the disease for a patient. DCMs in this context can already provide much less entangled voxelwise microstructure parameters such as the free water proportion (able to directly characterize edema) or intra / extra-axonal fractions at the voxel level. While these quantities can already be an improvement over DT, we have gone one step further [START_REF] Hédouin | [END_REF]] and provide a tract-based patient to population comparison framework for the characterization of changes within the patient microstructure at the compartment level.

Atlas-based patient to population study

Our framework is based on two parts: 1-an offline reference atlas creation providing all the necessary material about a control population for further patient study, 2the registration and patient analysis along reference fiber tracts on the atlas. Such an approach has the advantage of allowing a fast patient analysis requiring only an additional registration to perform the patient analysis.

Atlas construction and reference fiber tracts extraction

The first step of our framework, performed once for all patients analyzed is the creation of an atlas of reference tracts from a set of control subjects. Following precepts in previous studies [Goodlett et al. 2006], we have built an atlas from the diffusion images directly following the strategy highlighted in Section 4.2.2.2:

• Create an atlas from DT images using a method inspired from [Guimond et al. 2000] using the log-Euclidean framework [Arsigny et al. 2006b] and a DT registration method [Suarez et al. 2012] We have therefore collected data for this patient following the same protocol that was used for the control subjects in Section 4.2.2.2, i.e. dMRI with 30 directions on a b-value shell of 1000 s.mm -2 , image size of 128×128×55, voxel resolution of 2×2×2 mm 3 . We estimated DCMs from this data with the same parameters as for the controls, and compared them along fiber tracts of the left and right CSTs in the atlas. We report in Figure 4.5 the evaluation of differences along the left and right CSTs between the patient and controls for compartment AD. This figure shows differences inside the lesion region (highlighted in green on Fig. 4.5.c) along the tract for AD while these differences do not appear in the controlateral fibers (right CST). While other detections are visible here and there in fibers, their concentration is more prominent in that region, suggesting a lower value of AD for the patient. This highlights a change specific to the compartments along these fiber tracts, which is much more interpretable than global measures at the voxel level. We then report in Fig. 4.6 the same evaluation for compartment FA. Contrarily to AD, FA along the fibers does not vary significantly inside the lesion. This behavior is different from what is usually seen on voxelwise studies where FA and AD tend to vary together. Overall, these results are very interesting and suggest the feasibility of the application of our approach. Work remains however to be accomplished for a complete interpretation of those results since the data used is not ideal for DCM microstructure parameters estimation.

Combining relaxometry and diffusion for disease characterization

The final contribution of this chapter goes in a different direction. It highlights the potential of combining DCM and RCM for a better comprehension of diseases.

In this case, we studied the detection of Gadolinium (Gd) active lesions without using Gd. Active lesions in MS are a crucial marker of the disease activity, both for diagnosis [START_REF] Thompson | [END_REF]] and for the evaluation of a treatment. When present, they highlight that the disease is active, thus meaning that the current treatment is not working well enough to stop the disease from progressing. Monitoring active lesions is thus crucial to decide on the treatment adaptation for the patient. The classical way to assess active lesions is through the injection of a paramagnetic solution -the well known Gadolinium (Gd) -that highlights the blood brain barrier breakdown. While very practical, Gd however causes problems and is not indicated from some persons (allergy, kidney problems, etc.). Moreover recent studies [Gulani et al. 2017] have shown, without proving toxicity, that repeated injections of Gd tended to create Gd deposition in the brain. There is therefore a great interest to replace Gd injection by a less invasive method to assess active lesions.

On the other hand, we have seen in our research that both diffusion and relaxometry are sensitive to myelin and microstructural changes in the brain in very complementary ways. We have therefore investigated [Chatterjee et al. 2018b] a way to learn patterns (combining diffusion and relaxometry) of Gd active lesions in the brain to enable Gd lesion detection without injecting Gd.

Machine learning scheme for Gadolinium lesion detection

Our machine learning scheme for Gd lesions detection follows the scheme presented in Figure 4.7. This scheme is split into two parts. First we learn the features that characterize Gd lesions from a set of DCM and RCM images where the ground truth of Gd lesions have been manually delineated. Then, we apply this scheme to a new patient to detect within his lesions where active voxels are located. More precisely, the classification problem we have introduced starts from segmentations of T 2 hyperintense lesions (active or not), that can be obtained either manually or automatically. From these lesions and the ground truth, the training part chooses a set of non active voxels and active voxels randomly from the training set so that their numbers are balanced. From these voxels, feature vectors are constituted and used to train a non linear Support Vector Machine (SVM) classifier [START_REF] Cortes | [END_REF] best suited to classify those input voxels. Since this random selection of voxels may not be enough to detect correctly all active lesion voxels, this training process is repeated 100 times with randomly selected voxel feature vectors from the training set, leading to a set of 100 classifiers. In the testing 4.4. Relaxometry and diffusion for disease characterization 59 phase, we start from a segmentation of T 2 hyperintense lesions (active or not) and test each voxel feature vector against each of the 100 trained classifiers. If the majority of the classifiers assign the feature vector to the active lesion class, then the voxel is declared as active lesion. Note that probabilities of active lesions could also be given instead of a binary decision.

Quantitative MRI features

We have studied two types of quantitative MRI feature vectors, coming from either DCM images or RCM images. For DCM, we have used a simplified model with one free water compartment (diffusivity of 3×10 -3 mm 2 .s -1 ) and three directional compartments modeled as zeppelins [START_REF] Panagiotaki | Compartment models of the diffusion MR signal in brain white matter: a taxonomy and comparison[END_REF] with fixed radial diffusivity. We chose such a model as the dMRI acquisitions from which we estimate DCMs are clinical 30 directions / one b-value shell data (1000 s.mm -2 ). From these models (estimated using the method proposed in Chapter 2), we extracted voxelwise features more specific than the regular single tensor: free water weight, directional FA (i.e. weighted average of FAs over the anisotropic compartments), directional ADC and directional AD. These features are then stacked per voxel into a diffusion feature vector F D .

Similarly to diffusion, we constructed feature vectors from RCM images (estimated using Chapter 3), providing at each voxel a set of three fractions characterizing short, medium and large T 2 relaxation. Again, we stack those values into a vector F R at each voxel. From those two sets of vectors, we also tested a combined feature vector F RD = {F R , F D }.

Main results

I report here two main results of our experiments, more details being provided in [Chatterjee 2018, Chatterjee et al. 2018b]. The first experiment evaluated the interest of combining diffusion and relaxometry in this detection task by evaluating the overall accuracy, true positive rate and true negative rates in detecting the class of the voxels in the left out part of the training set at each stage (see Fig. We have then evaluated the capacity of our model learnt on the training set, with F RD as a feature vector, to correctly detect Gd positive parts of T 2 hyperintense lesions in a patient left completely out of the learning scheme. We applied the testing process to the patient's images and obtain the detections illustrated on three different axial slices in Fig. 4.9. This figure shows good prediction results of whether a voxel in the lesion is Gd positive or not. Quantitatively, the Dice score [START_REF] Dice | Measures of the amount of ecologic association between species[END_REF] for Gd positive lesion prediction is 0.64 and the one for Gd negative lesion prediction is 0.86. This indicates that, even if there is room for improvement, especially since the acquisitions used for this study are not ideal for microstructure parameters estimation, the detection rates obtained suggest a great potential for Gd lesions detection without the injection of Gd.

Conclusion and perspectives

Registration and processing of complex diffusion models

We have developed new algorithms and processing tools to accommodate the fact that diffusion models are a bit special in the field of computational anatomy. These are however crucial for defining new tools assessing the brain microstructure in var- ious diseases. While processing DT images has been well explored, no methods or only a few were available for the interpolation and processing (e.g. parameters extraction along fibers, tractography, etc.) of DCM images. Since this field has been prolific in terms of number of models, we have proposed a novel interpolation and averaging method, generic enough to be virtually applicable to any model encountered as long as two steps can be defined: a distance between two compartments, and a way to merge two compartments into one. While this method has been successful in our experiments, an investigation remains to be done on the nature of the interpolation that is desired by the user. As rightfully depicted by [Taquet et al. 2015], one may want (e.g. in the case of tractography) to preserve directions rather than microstructure parameters, one may otherwise want to preserve the signal simulated from the model (as we did), or more to preserve microstructure parameters.

The development of such extensions for our method only needs changes in the two aforementioned steps, but care must be taken depending on the application being aimed.

Additionally a direct and natural extension of this work, that we have started looking into in [START_REF] Commowick | [END_REF] (based on the registration framework in Appendix B), is the registration of DCM images. Once the interpolation is defined, only the similarity metric between models remains to accomplish this task. Again making such measures as agnostic as possible to the type of compartments will have to be at the center of such research. While such a registration is doable and attractive in terms of developments, studies will then be needed comparing methods registering different types of images (structural, DT or DCM images) to see the real benefit of complex model images registration and how to provide the best results in the best achievable timings as image resolution keeps on increasing with new scanner capabilities (e.g. the HCP diffusion data).

Finally, another topic of interest on this specific field of interest is tractography. While only discussed here, it is clear that this field requires great improvements [Maier-Hein et al. 2017, Jeurissen et al. 2019] to remove false positives, to which DCM based tractography can contribute. In addition, combining tractography with microstructure information, as started by [Daducci et al. 2015, Girard et al. 2017], coming not only from diffusion but also from other quantitative MRI modalities could be very valuable.

Combined analysis of quantitative MRI

Our experiments with the two PhDs of Sudhanya Chatterjee and Renaud Hédouin have demonstrated the great potential of using microstructure information 1-at the fiber level to detect changes of diffusion properties inside MS lesions ; 2-to detect, using a combination of diffusion and relaxometry imaging, lesions that are active without having to inject contrast agents such as Gd, for which the long term effects of their deposition in the brain are still unknown. These results, although preliminary, contribute to this objective. One direct improvement of the proposed methods will be to consider microstructure parameters as functions along the fibers and apply specific statistical schemes to them [START_REF] Goodlett | [END_REF]]. For the classification scheme, we are currently looking more in-depth at the influence of each feature on the results and towards providing for uncertain lesions, secondary decision schemes to further help the classification.

Larger studies will be needed to confirm this potential but it already directs us towards the definition of measures for evaluating the disease status, response to treatment, etc. Many works will be needed in this field, ranging from the definition of detection and characterization frameworks robust to all sorts of errors (image artifacts of various sources, registration to templates when lesions are present, etc.) to methods able to account for the temporal dimension as many neurological diseases are slowly evolving and require longitudinal studies for better detection. In the end, such studies may very well enable the early detection of subtle microstructure changes leading to a specific trend of a disease.

Chapter 5

General conclusions and perspectives

In this manuscript, I have covered the methodological developments we have accomplished towards the use of quantitative MRI in a clinical context, in order to provide more interpretable, reproducible measurements of the microstructure integrity in the brain. We have accomplished already a lot on various topics. We have first defined multi-compartment models for both relaxometry and diffusion, with a nice similarity between the two models. Exploiting this similarity, we have proposed maximum likelihood estimation methods for these models using the variable projection framework, providing relatively fast estimation of such complex models. Based on this estimation, we have proposed processing methods to exploit the models in disease studies. Finally, we have started working on applications of these models towards better disease understanding and the definition of patient specific measures of disease burden. These studies, although preliminary and on relatively small databases, showed very encouraging results: more specific and interpretable conclusions, ability to detect properties of MS lesions without contrast agents. For all these methods, due to the large variety of models, I have always put an emphasis on making them as agnostic as possible to the internal compartment definitions so that they can be applied to virtually any multi-compartment model. Also, I have always insisted on making all of the methods available open-source so that anyone from the field can access them and test our articles on their data. This lead to two repositories on Github: Anima1 and Anima scripts2 that contain all developments in this manuscript and more in the future.

Perspectives in diffusion and relaxometry Many future directions of research have already been mentioned in the previous chapters on dMRI and relaxometry. Going further on more general topics for these two modalities, I have forgotten so far about the acquisition side of things. The current clinical acquisition accepted nowadays for dMRI (30 gradient directions on one b-value shell) remains very limited for multi-compartment models. While we can push on making the best out of it (by imposing priors or model simplifications), one other way concerns the push of new acquisitions to the clinic, such as e.g. the CUSP acquisition from [Scherrer & Warfield 2012] that allows multiple b-values and more directions, especially with new scanner capabilities such as simultaneous multi-slice acquisitions. Going towards those however requires to convince clinicians of the actual importance of the models developed here. This is why we are currently working on applying these techniques to several diseases and introducing new acquisitions together.

On the side of relaxometry, the problem is even more pronounced. The methods and results I have presented are very encouraging towards applying these techniques as standard tools for neurodegenerative diseases evaluation. However, many acquisition developments are still needed to enable clinically compatible acquisitions for T 2 relaxometry. We have been collaborating a lot in the past few years with the Children's Hospital Boston on those aspects but work remains to be done to obtain artifact free, fast enough to acquire, sequences. Current acquisitions still require more than 10 minutes to obtain good quality relaxometry. Applying techniques such as compressed sensing or super-resolution algorithms e.g. [START_REF] Meurée | [END_REF]] could allow such clinical applicability.

As a final remark on the acquisition side, the methods that have been proposed for DCM estimation could very well be extended to benefit from new multidimensional dMRI [START_REF] Topgaard | Multidimensional diffusion MRI[END_REF] i.e. acquisitions with arbitrary diffusion gradient waveforms. This would simply require to change in Chapter 2 the b-value by a bmatrix without changing the core of the algorithm. We would then benefit from these new acquisitions very promising to get direct microstructure parameters such as axonal diameter [Drobnjak et al. 2010], thus leading to even more interpretable results on patients.

Finally on these perspectives centered around dMRI and relaxometry, there is now a need to combine those modalities beyond just detection tasks (as we did for Gd in Chapter 4). While this is very interesting in the first place, it will now be better to switch to purely descriptive, comprehensive analysis of patients with both modalities to derive more complete disease burden scores for patients, and also to build atlases of the normal population and its evolution.

General perspectives Going on some larger perspectives, all these topics have proven their value for disease studies, mainly centered on MS. One of my aims is now to apply these generic techniques to more diseases, each of which will surely bring their new challenges. Among them, dementia is a pathology that we have started working on where the brain connexions are changed both in terms of connexion density and microstructure properties in the brain. This will probably require more global disease burden studies including connectome analyses, with the challenge of integrating several microstructure properties. Additionally, one may mention traumatic brain injury or stroke, where additional constraints due to potentially large lesions in the brain will bring challenges on matching or following structures. However, it will probably help having a more in-depth evaluation of damage due to the lesions and hopefully in recovering brain function in the future.

Another topic of my research, that has been running for several years now, is pediatric research. We are now at the point where, thanks to advances in acquisitions, some if not all of the quantitative modalities could be applied to children or babies. It could bring great advances in the future to improve the robustness of the methods presented here so that we can apply multi-compartment techniques to neonates and children and better evaluate e.g changes due to prematurity. This will however require 1-more robust estimation especially to the smaller amount of data (directions or echoes) due to time constraints and to movement artifacts (e.g. [Mangin et al. 2002, Chang et al. 2005, Sairanen et al. 2018]); and 2-the development of longitudinal atlases either from cross-sectional or longitudinal data (research we have started to work on with the current PhD of Antoine Legouhy [START_REF] Legouhy | [END_REF]). Separated from studies on different diseases that I believe would bring clear advances, another point to achieve is longitudinal analysis. Many of the brain diseases are neurodegenerative i.e. the disease is slowly evolving inside the brain. With the arrival of new longitudinal databases such as HCP lifespan, it will now be possible to investigate the parallel evolution of patients and healthy controls on quantitative images, which could bring new information able to detect pathological evolution earlier. This could be achieved by applying what has been proposed for shape analysis or structural images analysis (e.g. parallel transport [Lorenzi & Pennec 2014, Cury et al. 2019]) to diffusion images in the first place and then other modalities as well.

While dMRI and relaxometry present nice properties to highlight microstructure changes, they are not the only quantitative modalities that may be obtained from MRI. Future studies could also benefit from other complementary images such as QSM providing maps of magnetic susceptibility, due among other factors to the presence of iron which is a crucial marker in MS, Parkinson's or Alzheimer's diseases. Arterial Spin Labeling is also a modality of choice for perfusion. In the future, having all these modalities in a joint evaluation framework, where the right modalities are picked for the disease under consideration, will be very valuable. The main problem or such a framework will be to generalize the current fiber bundle studies based on only one microstructure parameter to several parameters (as already mentioned for the connectome studies).

Finally, one limitation of all studies presented here is the relatively small number of patients or controls being studied. While interesting as a proof of concept for the method, large databases integrating all these quantitative modalities will be necessary to verify the applicability of the methods in real cases. The constitution of databases with one quantitative modality or two are currently under way (e.g. the HCP databases or UK biobank) but large databases with a nice set of quantitative images remain to be gathered. On these large databases, the generalization of atlasbased comparison will pose new problems due to the number of images: how to compute atlases in a reasonable time and expand them as images arrive? Should we move to multi-atlas analysis but then how to compare one image to so many? How to ensure registration between all these images is good enough so that average images are not too blurry (should we use Large Deformation Diffeomorphic Metric where I 3 is the 3×3 identity matrix and t i is a R 3 vector. In addition, its matrix logarithm is defined as a null matrix with a translation equal to that of A i (t i ). This transformation has proven to be very useful and is also the fastest to estimate. Originally, these three parameters have been interpreted as a sliding of the block along the grid of voxels of the second image. The easiest optimization method for a local block is thus a discrete grid search on the voxel grid of the floating image, which has the advantage of not requiring any interpolation (and is thus fast). Although enough to recover globally sub-voxel precision thanks to the iterations of the blockmatching algorithm, such an optimization may miss fine sub-voxel displacements.

Recent algorithms [Commowick et al. 2012b] have therefore studied optimization over the whole R 3 space, which may use a gradient-based algorithm if gradient of the similarity measure is available or gradient-free optimization such as the BOBYQA algorithm [Powell 2009].

B.2.2 Rigid transformation and beyond

Among linear transformations, virtually any can be searched for between a block and an image. One particularly interesting and that we have explored in several recent works [Commowick et al. 2012b, Commowick et al. 2012a] is the rigid transformation. It has several interesting properties that allow for its simple optimization and that will be useful to estimate global linear or non linear transformations. First of all, in homogeneous coordinates, A i is represented as follows:

A i = R i t i 0 1 (B.2)
where t i is again a translation (also integrating the center of rotation), and R i is a rotation matrix i.e. R i R T i = I 3 and ||R i || F = 1. Interestingly, Rodrigues' formulas allow for the explicit computation of its logarithm and exponential [Blanco 2010]. In particular, the matrix logarithm of a rigid matrix A is expressed as:

log(A) = w × l 0 0 (B.3)
where w × is the cross-product matrix of a vector of rotation angles w = (w 0 , w 1 , w 2 ) T and l an arbitrary R 3 vector. For the explicit logarithm and exponential formulas, refer to [Blanco 2010]. The parameterization of the transformation through the matrix exponential of log(A) is particularly interesting as it allows for 1-a clear and separate depiction of the six degrees of freedom (three scalars of l and the angles in w) of the transformation (instead of the rigid rotation matrix where parameters are entangled), and 2-a direct expression of the rigid transformation in its Lie group structure, useful to perform transformation extrapolation.

Coming back to block-matching, one can now search for rigid transformations between a block and the floating image instead of just looking for translations. This is particularly useful for registration of images with articulated structures such as the spine [Commowick et al. 2012a]. However, since the search space is now much larger (six variables instead of three), a brute force discrete search as in the translation case is much longer and, depending on the similarity measure optimized, a gradient-based or gradient-free optimization is more adapted.

As a final remark, we presented here the search over rigid transformations but any other transformation with well defined parameters covering its entire spectrum may be used (e.g. similarity transformation or more specific transformations as presented in Section 2.2).

B.3 Global transformation extrapolation

We now have from the previous step (actually the real block-matching step), a set of blocks B i defined on an image R and corresponding linear transformations Âi best matching them onto F . We now briefly detail how to use these local transformations to infer the update transformation δT . This procedure depends on the nature of local transformations and the nature of the global transformation sought.

B.3.1 Global linear transformations

Let us consider that we are looking for a global transformation that is linear. Two cases are still possible. First, if the Âi are translations, it is relatively easy to estimate any linear transformation best summarizing the local displacements obtained. Particularly, Chapter 8 of [Pennec 1996] explains very well how, through a linear least squares formulation, to estimate from a set of translations either a global translation, rigid or affine transformation and we refer the reader to these formulas for more details.

The second case arises when transformations Âi are more than just translations. In that case, we resort to the log-Euclidean framework for linear transformations [Arsigny et al. 2009] and thus formulate the least squares optimization directly on the matrix logarithms:

log(δT ) = arg max M i || log( Âi ) -M || 2 (B.4)
This directly leads to log(δT ) being the weighted average of the input logtransformations. For both cases, it is important to note that block-matching is not exempt from outliers which may degrade the obtained transformation. Numerous options may be applied to deal with this problem including weighting the individual Âi by the optimum similarity measure value they obtained, or performing robust estimation such as least trimmed squares or M-estimation [Rousseeuw & Leroy 1987].

B.3.2 Non linear transformations

When non linear transformations are sought, the extrapolation step is more difficult since many parameters come into play i.e. one 3D vector per voxel. In this section we will consider diffeomorphisms defined by their SVF as the final transformation being sought. While encompassing a reduced set of the diffeomorphisms that may be encountered (contrarily to LDDMM [Beg et al. 2005]), they have interesting properties that will be heavily used in this section, and that were well defined in the log-Euclidean frameworks for diffeomorphisms [Arsigny et al. 2006a] and polyaffine transformations [Arsigny et al. 2009].

In particular, we consider that δT is defined by its SVF δS i.e. a vector field whose exponential is δT : δT = exp(δS). The goal is to extrapolate the "log-vector" at each voxel of δS from the sparse set of optimal block transformations Âi located at their block centers x i . For this task, we will heavily use the matrix logarithm of linear transformations, which is explicitly defined for translations and rigid transformations.

As a side note to this class of transformations, it has been demonstrated that computing the SVF from a transformation T is computationally expensive [Arsigny et al. 2006a]. In addition, SVF having nice properties for statistics computation, it is desirable to always keep T implicit and instead compute the SVF S. Transposing the transformation composition in that space however requires to use the Baker-Campbell-Hausdorff (BCH) approximation [Bossa et al. 2007, Vercauteren et al. 2008].

B.3.2.1 Gaussian extrapolation

The first and simplest way to extrapolate a dense SVF is to use Gaussian extrapolation [START_REF] Commowick | Design and Use of Anatomical Atlases for Radiotherapy[END_REF], Garcia et al. 2010]. From the set of log( Âi ) and transformations locations x i , we first construct a sparse field C where each voxel corresponding to x i is affected by the displacement generated by log( Âi ): C(x i ) = log( Âi )x i . The Gaussian extrapolation then builds a dense SVF δS from C and a sparse field W of weights w i attributed to each pairing (for example the optimal value of the similarity measure): W (x i ) = w i . The extrapolated δS is then defined as:

δS = G σ * W C G σ * W (B.5)
where σ is the standard deviation of the Gaussian extrapolation kernel G σ . This extrapolation works perfectly in a region where enough input matchings are present, e.g. inside the brain. On the contrary, in regions far away from the blocks, this extrapolation is meaningless and may lead to artificially large deformations. To counter this effect, an additional post-processing is performed on the obtained SVF: when far enough from any matching (G σ * W below a certain threshold), δS is gradually set to identity (i.e. zero velocity field).

As for linear transformation computation, outliers in pairings need to be accounted for. A simple operation to do so is to compute at voxel locations x i the norm of the difference between δS(x i ) and C(x i ): r i = ||C(x i ) -δS(x i )||. From these, the mean displacement difference r over the whole image is computed as well as the variance σ 2 r :

not play the same role and reverting their use does not lead to the exact inverse of δT . Options are however available to ensure this, and this is why we detail three algorithms, going from no symmetry to more and more symmetry.

B.4.1 Asymmetric registration

The first, classical, registration built around the block-matching is the asymmetric one. It is described in Algorithm 1.

Algorithm 1 Asymmetric Block-Matching Registration Algorithm for l = 1...L, iterations, do 3:

Resample F with T

4:

Match R and F • T : δT ← block-match(R, F • T )

5:

Update T by composing it with δT 6:

If needed, regularize T (elastic-like)

From two images, a reference R and a floating image F , the algorithm seeks T by running a multi-resolution pyramid. At each step, the previously described components are put together to estimate update δT (or δS if the transformation computed is non linear) and compose it with the current estimate of T (BCH approximation for SVF). In this case, R and F clearly have an asymmetric role, blocks being always defined on R and only F being resampled.

B.4.2 Symmetric registration

In the previous algorithm, blocks are always defined on R while F is always the floating image. This second algorithm, presented in Algorithm 2, aims at symmetrizing this definition of blocks and ensuring that the obtained transformation when inverting F and R roles is the same up to an inverse, hereafter called inverse symmetry.

Algorithm 2 Symmetric Block-Matching Registration Algorithm Compute the transform update δT from δT F and δT R

7:

Update T by composing it with δT 8:

If needed, regularize T (elastic-like) Again, T and δT are replaced by S and δS when dealing with non linear transformations. In this algorithm, blocks are defined at each step both on R and on F and used to estimate two asymmetric updates: δT F and δT R (respectively δS F and δS R for non linear transformations). To account for these two updates and ensure inverse symmetry, the composition step is modified and preceded by an averaging of the two updates (for linear transformations using the matrix logarithm, and for SVF the log-Euclidean framework): Compute the half transform update δT from δT F and δT R

7:

Update T by composing it with δT 8:

If needed, regularize T (elastic-like)

As for direct symmetry in Algorithm 2, the composition step is modified to compute δT , respectively δS, from the asymmetric updates:

δT = exp 1 4 [log(δT R ) -log(δT F )] (B.11) δS = 1 4 (δS R -δS F ) (B.
12)

The only difference with direct symmetry is here 1/4 instead of 1/2 to account for the fact that we are looking for a transformation bringing the two images on a middle point where they match. Applying the final transformation T to F , is as simple as taking the square transformation (or multiply it by 2 in the "log-space"

Figure 1 . 1 :

 11 Figure 1.1: Illustration of the axon and myelin degeneration process (a) occurring in MS, leading to increased handicap on the Expanded Disability Status Scale (EDSS). The variability of this evolution is illustrated in (b,c).

Figure 2 . 1 :

 21 Figure 2.1: Illustration of several DW images with (a): no diffusion weighting, (b): diffusion gradient along the left-right axis, (c): diffusion gradient along the topdown axis. The red box illustrates the corpus callosum region, a region known for containing left-right fiber bundles. The blue box illustrates a part of the CST, a region known for containing top-down fiber bundles.

Figure 2 . 2 :

 22 Figure 2.2: Illustration of the different cell types of the brain: supporting (glial) cells and the neuron.

  (a) Left-Right (b) Right-Left (c) Block-matching (d) T1-weighted (e) Anterior-Posterior (f) Posterior-Anterior (g) Block-matching (h) T1-weighted

Figure 2 . 3 :

 23 Figure 2.3: Illustration of BM EPI distortion corrections on b 0 images acquired with opposite PEDs on one subject. First row: PEDs along the left-right axis, second row: PEDs along the anterior-posterior axis. (a-b, d-e) uncorrected EPI with opposite PEDs; (c,f) corresponding BM corrected images; (d,g) T 1 -weighted reference image.

Figure 2 . 4 :

 24 Figure 2.4: General illustration of the model estimation problem in dMRI.

Figure 2 . 6 :

 26 Figure 2.6: Illustration of the first elements of the SH basis. Images are courtesy of [Leemans 2010].

Figure 2 . 7 :

 27 Figure 2.7: Illustration of a DCM. First row: abstract model showing the combination of some isotropic compartment and several anisotropic compartments. Second row: illustration of the first row with a multi-tensors model.

Figure 2

 2 Figure 2.8: Illustration of the synthetic dataset used for ML estimation experiments. (a): Number of anisotropic compartments per voxel (dark blue: 0, cyan: 1, orange: 2, red: 3). (b): glyphs illustrating the DCM PDF.

Fig. 1 .Fig. 2 .Fig. 1 .Fig. 2 .Fig. 1 .Fig. 2 .Fig. 1 .Fig. 2 .Fig. 1 .Fig. 2 .Fig. 1 .Fig. 2 .Figure 2

 1212121212122 Fig. 1. Performance curves. The curves show the drop in MSE as a function of computation time (in sec). 1st row for the weights estimator, 2nd row for the tensors estimator. Columns represent areas with 0F, 1F, 2F, 3F from left to right.

Figure 2 .

 2 Figure 2.10: Visualization of the estimated multi-tensor models in vivo on an HCP subject. (a,b): estimation in the centrum semi-ovale, (c,d): estimation in the corpus callosum. Estimations are performed using method A1 (b,d) or method B (a,c).

Figure 3 . 1 :

 31 Figure 3.1: Illustration of relaxation times in MRI. T 1 : longitudinal relaxation time (a), T 2 and T * 2 : transverse relaxation times (b). Images are courtesy of [Ridgway 2010].

Chapter 3 .Figure 3 . 2 :

 332 Figure 3.2: Multiple T 2 components of a white matter voxel. (a): illustration (courtesy of [MacKay & Laule 2016]) on an electron microscopy image of the different components (red arrow: myelin, orange arrows: intra and extra-cellular matter). (b): illustration of a typical brain white matter voxel and its three T 2 components (short T 2 , medium T 2 , high T 2 ). (c): ratios between the different compartments.

CHAPITRE 3 .Figure 3 . 2 :

 332 Figure 3.2: Exemple de courbe réelle (en bleu) et de courbe estimée (en rouge) avec modèle exponentielle

Figure 3 . 4 :Figure 3 . 3 :

 3433 Figure 3.4: Exemple de courbe réelle (en bleu) et de courbe estimée (en rouge) avec modèle corrigé EP G

Figure 3 . 4 :

 34 Figure 3.4: Illustration of the numerical phantom developed for the numerical evaluation of our RCM model estimation framework.

FAE

Figure 3 . 5 :

 35 Figure 3.5: Relative mean square errors between ground truth values and estimated ones for estimated water fractions (left) and B 1 factor (right). Both graphs are a function of SNR (log-scale).

2 Figure 3 . 6 :

 236 Figure 3.6: Relative absolute differences between 32 and 7 echoes sequences for (a) short T 2 , (b) medium T 2 , (c) long T 2 water fractions. All errors are shown by region of interest as shown in Fig. 3.4.

Figure 3 . 7 :

 37 Figure 3.7: Illustration of the MS lesions study plan. Lesions (active: E+, and inactive parts: L-) are marked on the first time point and followed over 3 years.

Figure 3 .

 3 Figure 3.8: Illustration of the evolution of lesion water fraction values for an active lesion (at the first time point of the study) of a patient.

2 Figure 3 . 9 :

 239 Figure 3.9: Evolution of E+ and L-water fraction weights over the ten patients of the database and the three years. Stars indicate significant differences between E+ and L-at given time points (Mann-Whitney U test, p < 0.05).

Figure 6 . 5 :

 65 Figure 6.5: Test retest scans were performed for 4 healthy controls to study the repeatability of the proposed method. This figure shows the 15 regions which were marked on the healthy controls over which the repeatability was studied. The details of the acquired data are as follows: 3T MRI scanner; 2D multislice

Figure 3 .

 3 Figure 3.11: Water fractions repeatability over 15 different regions of interest on a test-retest experiment. (a,b,c): short, medium and high T 2 water fractions.

Figure 3 .

 3 Figure 3.12: Illustration of the estimated water fraction maps and medium T 2 mean relaxation time on an MS patient.

Figure 3 .

 3 Figure 3.13: Illustration of the medium T 2 compartment average relaxation time going through three different lesions, illustrating different lesion patterns.

Figure 4

 4

Figure 4 . 1 :

 41 Figure 4.1: Global interpolation scheme of DCM as a simplification problem.

Figure 4 . 2 :

 42 Figure 4.2: Illustration of DDI interpolation methods after three consecutive rotations around the same axis, compared to the original DDI. Top line: complete images. Middle line: zooms on the yellow region. Bottom line: zoom with compartments normalized to the same size to focus on the DDI compartments orientations.

Figure 4 . 3 :

 43 Figure 4.3: Example of a DDI atlas superimposed on the average b 0 image: axial view (first line) and coronal view (second line).

Figure 4 . 5 :

 45 Figure 4.5: Patient evaluation on the left and right CSTs for compartment AD. First line: left CST, second line: right CST. First column: patient AD. Second column: AD average across controls. Scalar bar below the two first columns: AD in mm 2 .s -1 . The last column corresponds to the p-value. Green volume: a lesion segmented manually on the T 2 -weighted image.

Figure 4 . 6 :

 46 Figure 4.6: Patient evaluation on the left and right CSTs for compartment FA. First line: left CST, second line: right CST. First column: patient FA. Second column: FA average across controls. Scalar bar below the two first columns: FA values. The last column corresponds to the p-value. Green volume: a lesion segmented manually on the T 2 -weighted image.

Figure 4 . 7 :

 47 Figure 4.7: Training and testing scheme for Gd lesions detection.

  4.7) by the same framework using either F D , F R or F RD . The results of this evaluation are presented in Figure4.8. This figure demonstrates the complementarity of the two modalities: diffusion alone allows a better true positive rate, on the contrary relaxometry alone allows a better true negative rate. When combining the two, both true positive and true negative rates reach values above each of the modalities alone. This translates to a better overall accuracy when combining the two modalities.

Figure 4 .

 4 Figure 4.8: Accuracy of predictions on the validation set depending on the features used for classification. From left to right: overall accuracy, true positive and true negative rates.

Figure 4 . 9 :

 49 Figure 4.9: Illustration of detection results on three different axial slices of an MS patient. Legend: blue: correctly predicted inactive lesions ; red: correctly predicted active lesions ; green: false positive voxels (i.e. detected as active when inactive).

1 :

 1 for p = 1...P , iteration on pyramid levels, do 2:

1 : 3 : 4 : 5 :

 1345 for p = 1...P , iteration on pyramid levels, do 2:for l = 1...L, iterations, do Resample F with T and R with T -1 Match R and F • T : δT F ← block-match(R, F • T ) Match F and R • T -1 : δT R ← block-match(F, R • T -1 ) 6:

3 Algorithm 3 2 : 3 : 4 : 5 :

 332345 Kissing symmetric registrationIn direct symmetry, the images roles are not completely symmetric. In fact, one image in each way (from F • T to R and from R • T -1 to F ) is never resampled: the transformation is applied only to one image at a time. Kissing symmetry instead seeks the transformation T so that R • T -1 and F • T match. T now encodes the half transformation between the images: this amounts to looking for an intermediate position in between the two images by moving both of them towards each other, thereby fully symmetrizing their roles. This registration is presented in Algorithm 3. Kissing Symmetric Block-Matching Registration Algorithm 1: for p = 1...P , iteration on pyramid levels, do for l = 1...L, iterations, do Resample F with T and R with T -1 Match R • T -1 and F • T : δT F ← block-match(R • T -1 , F • T ) Match F • T and R • T -1 : δT R ← block-match(F • T, R • T -1 ) 6:

  Introduction Contents 1.1 Foreword . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2 Quantitative images for disease understanding . . . . . . . .

	Chapter 1

2 1.2.1 An example: multiple sclerosis and MRI . . . . . . . . . . . . 2 1.3 Contributions summary . . . . . . . . . . . . . . . . . . . . . . 4 1.3.1 Diffusion imaging for white matter microstructure imaging . 4 1.3.2 Getting insights on myelin degeneracy: relaxometry . . . . . 4 1.3.3 Quantitative image processing for disease study . . . . . . . . 4

  Artifacts correction in diffusion weighted imaging . . . . . . 10 2.2.1 EPI and distortion artifacts . . . . . . . . . . . . . . . . . . . 10 2.2.2 Distortion model . . . . . . . . . . . . . . . . . . . . . . . . . 11 2.2.3 A block-matching strategy for distortion correction . . . . . . 11 2.2.4 Main results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.3 Modeling water diffusion . . . . . . . . . . . . . . . . . . . . . 15 2.3.1 Diffusion tensor . . . . . . . . . . . . . . . . . . . . . . . . . . 16 2.3.2 Orthogonal bases . . . . . . . . . . . . . . . . . . . . . . . . . 18 2.3.3 Diffusion compartment models . . . . . . . . . . . . . . . . . 19 2.4 Estimation of diffusion compartment models . . . . . . . . . 21 2.4.1 Estimation problem formulation . . . . . . . . . . . . . . . . 21 2.4.2 Variable projection solution . . . . . . . . . . . . . . . . . . . 22 2.4.3 Main results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 2.5 Conclusion and perspectives . . . . . . . . . . . . . . . . . . . 26 2.5.1 Artifacts correction in diffusion weighted imaging . . . . . . . 26 2.5.2 Diffusion models estimation . . . . . . . . . . . . . . . . . . . 27

	Chapter 2
	White matter microstructure from
	diffusion imaging
	Contents
	2.1 Diffusion imaging and white matter microstructure . . . . .	8
	2.2

  Chapter 2. White matter microstructure from diffusion imaging 2.1 Diffusion imaging and white matter microstructure

	, Olivier
	Commowick, Elise Bannier, Benoit Scherrer, Maxime Taquet, Simon K
	Warfield and Christian Barillot. Block-Matching Distortion Correction of
	Echo-Planar Images With Opposite Phase Encoding Directions. IEEE Trans-
	actions on Medical Imaging, 36(5):1106-1115, 2017
	• Diffusion model estimation: Aymeric Stamm, Olivier Commowick, Simon K. Warfield and Simone Vantini. Comprehensive Maximum Likelihood Estima-

tion of Diffusion Compartment Models Towards Reliable Mapping of Brain Microstructure. In MICCAI proceedings, pages 622 -630, 2016. Diffusion MRI (dMRI) [Le

  .1. These results have highlighted that BM improves the correction over Voss et al. on all subjects. Between BM and TOPUP, the best score depends on the subject with a similar average, highlighting similar performance.

		Correlation Untouched Voss	BM TOPUP
		Subject 1	0.842	0.901 0.916	0.927
		Subject 2	0.818	0.904 0.918	0.937
		Subject 3	0.812	0.875 0.894	0.859
		Subject 4	0.886	0.923 0.939	0.954
		Subject 5	0.872	0.913 0.921	0.898
		Mean	0.852	0.903 0.918	0.915
	Table 2.1:	Correlation results between left-right/right-left and anterior-
	posterior/posterior-anterior images.

  Contents 3.1 Quantitative relaxation times from MRI . . . . . . . . . . . 29 3.2 Modeling multiple T 2 compartments . . . . . . . . . . . . . . 32 3.3 Robust compartment models estimation . . . . . . . . . . . 34 3.3.1 Non-negative compartment weights estimation . . . . . . . . This work has been conducted mainly with one PhD student, one post-doc and one intern I co-supervised: Sudhanya Chatterjee, Fang Cao and Lucas Soustelle.Several papers resulted from this work, but particularly two main ones discussed in the following sections:

	• Robust Relaxometry Compartment Model (RCM) estimation and use: Sud-hanya Chatterjee, Olivier Commowick, Onur Afacan, Benoit Combes, Simon
	Warfield, Christian Barillot. A three year follow-up study of gadolinium en-
	hanced and non-enhanced regions in multiple sclerosis lesions using a multi-
	compartment T2 relaxometry model. Preprint on bioRxiv. 2018
	• Variable projection for RCM relaxometry: Sudhanya Chatterjee, Olivier Commowick, Onur Afacan, Simon K Warfield, Christian Barillot. Multi-
	Compartment Model of Brain Tissues from T2 Relaxometry MRI Using
	Gamma Distribution. In ISBI proceedings, pages 141 -144, 2018.

34 3.3.2 Main results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 3.4 Towards RCM parameters estimation . . . . . . . . . . . . . 41 3.4.1 Gamma relaxometry compartment model . . . . . . . . . . . 41 3.4.2 Maximum likelihood estimation framework . . . . . . . . . . 41 3.4.3 Main results . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42 3.5 Conclusion and perspectives . . . . . . . . . . . . . . . . . . . 45

This chapter explores our research around quantitative relaxometry sequences.

  .7.

	same regions are observed over 3 years
	M00 M03 M06 M09 M12	M18	M24	M36
		Illustration of lesion ROIs studied
	regions marked		appears on T2-w but not on Gd (L-)
	at M00
			appears only on
			Gd post contrast
			(E+)
		MS lesion on T2-w MRI	Regions studied

•

  DCM interpolation: Renaud Hédouin, Olivier Commowick, Elise Bannier, Aymeric Stamm and Christian Barillot. Interpolation and Averaging of Multi-

	Compartment Model Images. In MICCAI proceedings, pages 354 -362, 2015.

  ). & exhibition, Honolulu, Hawaii, USA, United States, April 2017. (Cited on page 34.) [Chatterjee et al. 2018a] Sudhanya Chatterjee, Olivier Commowick, Onur Afacan, Benoit Combes, Simon K Warfield and Christian Barillot. A three year follow-up study of gadolinium enhanced and non-enhanced regions in multiple sclerosis lesions using a multi-compartment T2 relaxometry model. bioRxiv, 2018. (Cited on pages 36, 40 and 43.) [Chatterjee et al. 2018b] Sudhanya Chatterjee, Olivier Commowick, Onur Afacan, Simon K Warfield and Christian Barillot. Identification of Gadolinium contrast enhanced regions in MS lesions using brain tissue microstructure information obtained from diffusion and T2 relaxometry MRI. In 21st International Conference on Medical Image Computing and Computer Assisted Intervention -MICCAI, volume 11072 of LNCS, pages 63-71, September 2018. (Cited on pages 58 and 60.) [Chatterjee et al. 2018c] Sudhanya Chatterjee, Olivier Commowick, Onur Afacan, Simon K Warfield and Christian Barillot. Multi-Compartment Model of Brain Tissues from T2 Relaxometry MRI Using Gamma Distribution. In ISBI -IEEE International Symposium on Biomedical Imaging, pages 141-144, 2018. (Cited on pages 41 and 42.) [Chatterjee 2018] Sudhanya Chatterjee. Gaining insights into brain tissues using multi-compartment T2 relaxometry models. PhD thesis, Universite Rennes 1, December 2018. (Cited on pages 36, 46 and 60.) [Chung et al. 2011] J.Y. Chung, M.H. In, S.H. Oh, M. Zaitsev, O. Speck and Z.H. Cho. An improved PSF mapping method for EPI distortion correction in human brain at ultra high field (7T). Magnetic Resonance Materials in Physics, Biology and Medicine, vol. 24, no. 3, pages 179-190, 2011. (Cited on page 10.) [Combès et al. 2016] Benoit Combès, Anne Kerbrat, Olivier Commowick and Christian Barillot. A Bayesian Model to Assess T2 Values and Their Changes Over Time in Quantitative MRI. In 19th International Conference on Medical Image Computing and Computer Assisted Intervention (MICCAI), pages 570 -578, Athens, Greece, October 2016. (Cited on page 31.) [Commowick et al. 2012a] Olivier Commowick, Nicolas Wiest-Daesslé and Sylvain Prima. Automated diffeomorphic registration of anatomical structures with rigid parts: application to dynamic cervical MRI. In 15th International Conference on Medical Image Computing and Computer Assisted Intervention, volume 15 of LNCS, pages 163-70, Nice, France, October 2012. Springer. (Cited on pages 74, 75 and 78.) [Commowick et al. 2012b] Olivier Commowick, Nicolas Wiest-Daesslé and Sylvain Prima. Block-Matching Strategies for Rigid Registration of Multimodal Medical Images. In 9th IEEE International Symposium on Biomedical Imaging

http://www.crl.med.harvard.edu

https://doi.org/10.17605/OSF.IO/YP4QG

http://github.com/Inria-Visages/Anima-Public, RRID:SCR_017017

http://github.com/Inria-Visages/Anima-Scripts-Public, RRID:SCR_017072

• Apply the obtained transformations to DCM images with our interpolation

• Compute an average DCM image.

From this atlas, we then extract the fiber tracts of interest using a full brain tractography adapted to work with DCM images. To do so, we adapted a classical deterministic tensor tractography algorithm [Mori et al. 1999] to handle DCM images in a generic manner (i.e. independent of the underlying compartment model). The algorithm is a seed-based one, and its main adaptation is to take the direction for the next progression step as the one from the closest DCM compartment to the previous direction. This allows especially for a better handling of crossing fibers in the brain. This tractography is, again, made once and for all on the atlas and will be used as a basis to extract compartment-based microstructure information for all control subjects and for the patient.

Patient fiber tracts analysis

Once the atlas of control subjects is built, the second step of the framework consists in registering the patient DCM image to be evaluated on it. This image is registered in the same way as for the atlas construction, to avoid bias in the results. With this registration performed, the fiber tracts of interest are used to extract microstructure properties along the fibers (as shown in Figure 4.4) both in the patient and the control subjects. Statistical tests can then be performed to characterize the difference between the patient and controls. This last part is performed as proposed in [START_REF] Commowick | [END_REF] where we test if the parameter of the patient is significantly different from the controls taken as a normal distribution.

Main results

As a first proof of concept of this new framework, we have evaluated in [START_REF] Hédouin | [END_REF]] a patient suffering from MS, with one lesion along the left CST.

Appendix A

Variable projection for maximum likelihood model estimation We present in this appendix the general derivation of the variable projection estimation framework for maximum likelihood estimation of models expressed as a linear combination of independent compartments. This derivation closely follows the work from [Golub & Pereyra 1973] and is used for both diffusion DCM estimation in Chapter 2 and RCM estimation in Chapter 3.

A.1 Problem formulation and maximum likelihood estimation

We consider the case of general estimation problems where a set of signals y = {y 1 , . . . , y N } T is acquired and modeled by a multiple compartment equation, i.e. a linear weighted sum of individual models each representing a subpart of the signal formation. In practice, we consider the following model:

where x is e.g. the q-vector in Chapter 2 or the echo time e in Chapter 3, S 0 > 0 is a constant baseline signal, w = {w 1 , . . . , w C } T is a set of weights with the property that w j ≥ 0 and j w j = 1, and f j is a, usually non linear, function parameterized by the set of parameters θ j = {θ j,1 , . . . , θ j,k , . . . , θ j,K j }. We now consider the maximum likelihood estimation problem of the function f knowing signals y. We first assume that y is perturbed by some noise, i.e.

where ε follows a noise distribution. In this work, we will consider only white Gaussian noise, i.e. ε ∼ N (0, σ 2 ), although other noise types such as Rician noise could be considered. From this noise assumption, we define the Gaussian log-likelihood function for the problem:

where τ 2 = 1/σ 2 . The objective in maximum likelihood approaches is then to maximize this equation with respect to its parameters w, S 0 and θ. Without loss of generality and to simplify constraints setting, we first reparameterize f to fuse w j and S 0 values into a single set of parameters α = {α 1 , . . . , α C } T where α j = S 0 w j . By construction, α j ∈ R + and values of S 0 and w j can be recovered from α as S 0 = j α j and w j = α j /S 0 . f is thus now written as:

Written in matrix form, Eq. (A.3) becomes the following:

where

A.2 Projecting linear variables of the system

From this point on, several variables can be identified as linearly separable in the optimization of the log-likelihood in Eq. (A.5). First of all, if we study its derivative with respect to τ 2 and equate it to zero, we find the following expression for the optimal τ 2 as a function of other parameters:

Putting back this analytical solution inside Eq. (A.5), and after simplification, we then obtain the following equivalent formulation of the log-likelihood:

Analyzing this equation, it appears that maximizing this likelihood is equivalent to minimizing the following least squares problem:

If we look closely at this least squares system, some variables are again linear in the system: α, while some are non linear: θ. For such cases, [Golub & Pereyra 1973] proposed the variable projection technique. It consists in expressing the system only in terms of the non linear variables of the system, the linear ones being directly A.3. Differentiation against non linear parameters 71 determined by their analytical solutions in terms of the non linear parameters. In practice, for the current system, we first determine the optimal α as the classical solution of linear least squares:

where F + denotes the pseudo-inverse of F i.e. F + = F T F -1 F T . Putting back this solution inside the previous non linear least squares system, we get the following variable projection least squares system:

where I N denotes the N × N identity matrix. From this system, the estimation process is thus as follows:

• Estimate θ by solving the variable projected system in Eq. (A.10)

• Deduce from θ the value of α using Eq. (A.9)

• Finally compute τ 2 using Eq. (A.6)

A.3 Differentiation against non linear parameters

The only step that remains is the optimization against the non linear parameters.

For some estimation problems, computing the derivatives of the cost function in Eq. (A.10) is too complicated or even not analytically feasible. In those cases, Eq. (A.10) can be optimized using a derivative free algorithm such as the BOBYQA algorithm [Powell 2009].

For other problems where the derivative of the cost function may be computed, any derivative based algorithm may be used such as the BFGS optimizer [Byrd et al. 1995] or, as the residuals can be individually separated in Eq (A.10), the Levenberg-Marquardt optimizer [Levenberg 1944]. Let us rewrite E in a way suitable for the Levenberg-Marquardt optimizer. We define r as the residuals vector for each measurement:

We may redefine E as the following function of r: E(θ = r(θ), r(θ) , where ., .

denotes the dot product. The derivative of each residual of the orthogonal projection of the system is given in [Golub & Pereyra 1973] as the following general formula:

where ∇r is an N × N C matrix representing the derivatives of r against the N C parameters inside θ. By some abuse of notation for simplification D F represents, for each of the N C parameters inside θ, the component-wise derivative of F against this parameter. With D F computed, the general cost function derivative can then be expressed as:

An interesting observation about those two last equations is that to compute the derivative of any estimation problem, with any model, only the derivatives of this model against the parameters (the D F matrices) need to be computed for the variable projection problem to be optimized. This fact makes the variable projection technique generic and this is heavily used for DCM estimation in Chapter 2 or RCM estimation in Chapter 3.

Appendix B

A versatile framework for images registration We present in this appendix the generic registration algorithm used in several chapters (for different applications and image types). This algorithm relies on the block-matching algorithm introduced by [Ourselin et al. 2000] for rigid registration. It is in fact very versatile, usually requiring only the similarity measure between blocks (and a re-orientation scheme for images of oriented models such as DT or DCM images) to be defined to perform the registration. We explore in the first sections the individual elements of the registration and provide the general algorithms in Section B.4.

Contents

B.1 Block-matching for medical images registration

We define the registration algorithm as the one that seeks a transformation T so that a floating image resampled by T , F • T , matches as much as possible a reference image R: R ≈ F • T . With that goal in mind, the core part of the registration algorithm is based on an iterative framework which iterates three main steps:

• Define a set of blocks, i.e. a subset of voxels -often a cube around a given point in the image, on the reference image R. Each block is defined by its center x i : B(x i ) ≡ B i

• Match these blocks from image R to image F , i.e. find the best local transformation A i such that a similarity measure S(R(x), F • A(x)) over the block B i is optimal

• From the set of transformations obtained, compute an update global transformation δT then used to update the current global transformation T

This algorithm is the block-matching core that is central to all variants that are used in this document. This core is included in an iterative algorithm (see Section B.4), itself often included in a multi-resolution (pyramidal) framework to first get back large displacements from coarse resolutions and then smaller / finer displacements from finer resolutions, in a robust and faster manner. This core part has been used in many algorithms ranging from linear [Ourselin et al. 2000, Commowick et al. 2012b] to non linear, diffeomorphic registration [Commowick et al. 2012a], with various applications such as tensor registration [Suarez et al. 2012]. It requires mainly the definition of a similarity metric S between blocks to be used on a specific modality. On scalar images [Malandain et al. 2004], given the relatively small size of a block (thus containing two or few different tissues), a linear relationship between intensities is usually enough and a square correlation coefficient is used to match images. Apart from the similarity, several points need to be defined at the general level to have an overview of the algorithm:

• The nature of transformations A and T

• How to go from a set of transformations A i and blocks B i to a global (linear or non linear) update transformation δT ?

• Ensuring (or not) a symmetric transformation

B.2 Which local transformations between blocks?

Local transformations between a block in R and F are generally assumed to be linear for two reasons: 1-the block is usually sufficiently small for this assumption to be true, and 2-the search space over which the transformation has to be estimated needs to remain small enough for computation time reasons. A transformation A i in R N is therefore represented as a (N + 1) × (N + 1) matrix. In the following, let us consider without loss of generality that N = 3.

B.2.1 Local translation

The most common transformation used in block-matching is a local translation, i.e. a move of the block in the three directions. In other words, A i is represented as:

B.2. Which local transformations between blocks

We can reject pairings from C for which the residual r i > r + ασ r and recompute from it an outlier free δS.

B.3.2.2 M-Smoothing extrapolation

Although simple, the Gaussian extrapolation does not completely incorporate outlier rejection in its process in the sense that one would need to iterate over the rejection process to do so.

We have therefore introduced in [Commowick et al. 2012a] an extrapolation approach similar to the M-smoothing filter proposed by [START_REF] Mrazek | [END_REF]]. This approach looks for the best logtransformation log(S k ) at each voxel of δS by minimizing the following criterion:

where ρ is robust error norm (typically linked to an M-estimator, here the Welsch function), V k is a neighborhood around voxel i (note that the sum over i ∈ V k only considers voxels where a transformation A i was estimated) and d is a spatial error norm. This criterion can be minimized using gradient descent which for a particular adaptive, data-dependent step size leads to the following update formula for each log S k :

where ρ acts as a tonal kernel, which for the Welsch function ρ is written as ρ (a 2 ) = exp -a 2 /2λ 2 , and d acts as a spatial kernel, here a Gaussian kernel: d(a 2 ) = exp -a 2 /2σ 2 . The gradient descent is initialized with ρ (a 2 ) = 1. These two kernels account simultaneously for the spatial proximity of A i and its local agreement with other local transformations. From the obtained log S k , we finally obtain δS at each voxel by δS(x k ) = log S k x k .

B.4 Asymmetric or symmetric registration

From the previous sections, we now have all the necessary elements (apart from a few such as being able to resample images or the similarity measure which is not the topic here) to perform the registration of two images. The final step is to combine all of these into an algorithm. At this stage, it is crucial to note that the blockmatching core algorithm is intrinsically an asymmetric one: images R and F do B.4. Asymmetric or symmetric registration