Renaud H Édouin
email: renaud.hedouin@gmail.com

Member, IEEE Christian Barillot
email: christian.barillot@irisa.fr

Olivier Commowick
email: olivier.commowick@inria.fr

Interpolation and averaging of diffusion MRI multi-compartment models

Keywords: Multi-compartment models, Interpolation, Diffusion MRI, Atlases

Multi-compartment models (MCM) are increasingly used to characterize the brain white matter microstructure from diffusion-weighted imaging (DWI). Their use in clinical studies is however limited by the inability to resample an MCM image towards a common reference frame, or to construct atlases from such brain microstructure models. We propose to solve this problem by first identifying that these two tasks amount to the same problem. We propose to tackle it by viewing it as a simplification problem, solved thanks to spectral clustering and the definition of semi-metrics between several usual compartments encountered in the MCM literature. This generic framework is evaluated for two models: the multi-tensor model where individual fibers are modeled as individual tensors and the diffusion direction imaging (DDI) model that differentiates intra-and extra-axonal components of each fiber. Results on simulated data, simulated transformations and real data show the ability of our method to well interpolate MCM images of these types. We finally present as an application an MCM template of normal controls constructed using our approach.

I. INTRODUCTION

Diffusion weighted imaging (DWI) is a specific type of MRI acquisition that studies the constrained random diffusion of water molecules within different tissues. It allows, through several acquisitions, to model the brain microstructure, i.e. white or gray matter internal cellular organization, which is significantly smaller than the voxel resolution. The most widespread model, especially in clinical routine, is named diffusion tensor imaging (DTI) [START_REF] Basser | Microstructural and physiological features of tissues elucidated by quantitative diffusion-tensor MRI[END_REF]. It represents the water diffusion within each voxel as a tensor, i.e. a zero-mean 3D Gaussian distribution. The DTI, simple to estimate and used in clinical studies for many diseases, also holds limitations. The diffusion in complex areas, such as in the white matter (WM) crossing fibers, cannot be represented by a single tensor. To overcome this issue, multi-compartment models (MCM) [START_REF] Ferizi | A ranking of diffusion MRI compartment models with in vivo human brain data[END_REF] (a class of diffusion models) consider separately groups of tissues with similar diffusion properties (e.g. mixture of main directions, diffusion pattern, etc.) and represent each of them by an individual compartment. For example, an anisotropic compartment models a fascicle, i.e. a bundle of axons with a similar general orientation.

MCM provide a greater specificity and sensitivity to microstructural properties than conventional DTI. Thus, it provides valuable information for a deeper understanding of the healthy brain. Moreover, it can be used in a clinical setting for providing new imaging biomarkers of inflammation, neuronal or axonal degeneration. Promising studies have shown that MCMs are efficient to distinguish between tumor types [START_REF] Panagiotaki | Noninvasive quantification of solid tumor microstructure using VERDICT MRI[END_REF], or detect specific tissue subtypes in multiple sclerosis [START_REF] Lakhani | Advanced multicompartment diffusion MRI models and their application in multiple sclerosis[END_REF]. A more general use of MCM in the clinic would greatly help for diagnosis or adapting treatment, in particular for diseases impacting WM such as neurodegenerative or neuroinflammatory diseases.

Drawing statistical conclusions from multi-subject DWI datasets is made difficult by the challenge of projecting DWI images into a common reference space, since they often do not share the same gradient directions or number of b-values. In addition, the statistical comparisons will finally be made on properties extracted from the MCMs. Doing the registrations or transformations application directly on the MCMs is a more direct solution. This requires some critical steps to be defined such as the application of a transformation to MCM images, the creation of MCM atlases, i.e. a common reference frame for comparison of patients and healthy controls, etc. Doing so requires however to define the interpolation task, i.e. finding the model in between voxels locations, and the averaging of several voxels into a single one, e.g. when merging several registered MCM images into the final atlas. These common steps for anatomical images are however not trivial for such complex diffusion models.

To date, only few approaches have addressed these issues for complex diffusion images including MCM. Among them, Barmpoutis et al. [START_REF] Barmpoutis | Registration of high angular resolution diffusion MRI images using 4th order tensors[END_REF] or Geng et al. [START_REF] Geng | Diffusion MRI registration using orientation distribution functions[END_REF] introduced registration methods specifically tuned for orientation distribution function (ODF) on the sphere. Goh et al. [START_REF] Goh | A nonparametric Riemannian framework for processing high angular resolution diffusion images and its applications to ODF-based morphometry[END_REF] introduced an interpolation method for ODF in a spherical harmonics basis as a Riemannian average. This is made possible by the fact that the spherical harmonics basis used is an orthogonal basis that allows the definition of a Riemannian manifold and thus classical operations such as the Frechet mean. However, this does not apply to MCMs. The closest approach to our work was proposed by Taquet et al. [START_REF] Taquet | Registration and analysis of white matter group differences with a multi-fiber model[END_REF] who introduced an interpolation approach seen as a simplification problem. Indeed, each voxel of the interpolated image is simply a weighted sum (8 in the case of a trilinear interpolation) of the neighbourhood voxels of the original image, each one being a weighted sum of compartments in the case of a MCM. In their method, a large number of weighted compartments from a set of voxels is summarized into a smaller set of compartments. However, they assume that a single compartment follows a simple exponential model, which is not the case for all MCM.

We therefore introduce a new interpolation and averaging method for MCM data with the goal of having a scheme as agnostic as possible to the particular combination of compartments used in an application. To do so, we first introduce interpolation and averaging as a simplification problem, after presenting in Section III-A these two problems and how they are the same problem. Our method relies on the general problem of the clustering from a set of input compartments into a predefined number of output compartments. Each cluster is then simplified into an interpolated compartment whose weight is computed from initial input weights. The combination of the interpolated compartments provides the output MCM. Unlike Taquet et al., this method is very generic as it relies only on the definition of a similarity measure between compartments and of a weighted averaging scheme for compartments. It can therefore be applied to any type of MCM compartments as long as those two components are defined.

After some background on MCM and the models used in our experiments in Section II, we first present the MCM simplification problem independently of the particular implementation of the MCM model in Section III. We then describe in Section IV the implementations of this global scheme for isotropic compartments and for two types of anisotropic compartments, tensors and diffusion directions imaging (DDI) [START_REF] Stamm | A new multi-fiber model for low angular resolution diffusion MRI[END_REF]. To this end, we define one metric for the tensor compartments and four possible metrics for the more complex DDI compartments. We demonstrate qualitatively and quantitatively the interest of our method on simulated and in vivo data. We finally apply this framework to compute an atlas of DDI models (Section V) that clearly highlights a better averaging of fiber crossing regions.

II. BACKGROUND ON MULTI-COMPARTMENT MODELS

A. Diffusion multi-compartment model formulation

A diffusion multi-compartment model (MCM) provides information, through a weighted sum of individual compartments, on the different elements of the brain microstructure at a given voxel in the brain. Each compartment therefore represents specific diffusion properties in a specific tissue (cerebrospinal fluid (CSF), grey matter or axons within a specific direction). Mathematically, a MCM probability density function (PDF) P(x) is expressed as a sum of isotropic and anisotropic compartments:

P(x) = M i=1 w iso,i p iso,i (x) + N j=1 w j p j (x) (1)
where p iso,i and p j are the PDFs of respectively the i -th isotropic compartment, and the j-th anisotropic compartment of the model. w iso,i and w j are the compartment weights of the model and sum up to 1.

Although this equation could be written in a simpler manner, we prefer to explicitly separate isotropic and anisotropic compartments as they represent totally different tissue types. More precisely, an isotropic compartment represents isotropic water diffusion inside the brain. This concerns water diffusion inside glial cells or neuron cell bodies, as well as free water diffusion within the brain. Various compartment models exist for these isotropic compartments [START_REF] Panagiotaki | Compartment models of the diffusion MR signal in brain white matter: a taxonomy and comparison[END_REF], [START_REF] Stanisz | An analytical model of restricted diffusion in bovine optic nerve[END_REF]. An anisotropic compartment models a fascicle, i.e. a bundle of axons with a similar general orientation. The behavior of water within such an environment is complex. A large number of models has therefore been proposed to describe these compartments. They sometimes distinguish between intra-and extra-axonal diffusion (e.g. NODDI [START_REF] Zhang | NODDI: practical in vivo neurite orientation dispersion and density imaging of the human brain[END_REF], DDI [START_REF] Stamm | A new multi-fiber model for low angular resolution diffusion MRI[END_REF]), or may be simpler (e.g. tensors [START_REF] Taquet | Registration and analysis of white matter group differences with a multi-fiber model[END_REF]).

Given the large variety of compartment models, both for isotropic and anisotropic ones, we have chosen here to present a versatile method that can adapt virtually to any of the compartment models. This general method is described in Section III-A. Since we need to evaluate this method on real MCM datasets though, we have chosen to evaluate it on two models described in the next sub-sections: the multi-tensors model (MTM) and the DDI model. Again these are models inspired from the literature and other options are possible both for the choice of isotropic compartments and anisotropic ones.

B. Multi-tensors model

The MTM, popularized for registration by Taquet et al. [START_REF] Taquet | Registration and analysis of white matter group differences with a multi-fiber model[END_REF], is a straightforward extension of the well known single tensor model. In our implementation, we have chosen to use a model composed of several elements. First we use two isotropic compartments: one to model free water diffusion, and one to model restricted isotropic water diffusion that happens in glial cells, neuron cell bodies [START_REF] Panagiotaki | Compartment models of the diffusion MR signal in brain white matter: a taxonomy and comparison[END_REF], [START_REF] Stanisz | An analytical model of restricted diffusion in bovine optic nerve[END_REF], and damaged tissues [START_REF] Wang | Quantification of increased cellularity during inflammatory demyelination[END_REF]. Both compartments are modeled with isotropic zero-mean 3D Gaussian PDFs, but with different diffusivities: d free = 3 × 10 -3 mm 2 .s -1 for free water [START_REF] Pasternak | Free water elimination and mapping from diffusion MRI[END_REF]- [START_REF] Harris | Pressure and temperature dependence of the self diffusion coefficient of water and oxygen-18 water[END_REF], and d free = 1 × 10 -3 mm 2 .s -1 for restricted water diffusion.

In addition to these isotropic compartments, anisotropic compartments are modeled as individual tensors with different main directions, i.e. zero-mean 3D Gaussian PDFs parametrized by their covariance matrices. In our experiments, we used three anisotropic tensors but this number can be smaller (or larger even though it is commonly admitted that at most three fibers directions cross in a single voxel [START_REF] Schilling | Can increased spatial resolution solve the crossing fiber problem for diffusion MRI?[END_REF]) depending on the MTM estimation scheme or if model selection was applied.

C. Diffusion directions imaging model

The second MCM we are running experiments on is the diffusion directions imaging (DDI) model. Similarly to the MTM, we model diffusion in cell bodies and free water by two isotropic compartments, each modeled by an isotropic Gaussian PDF with different diffusivities: d free = 3×10 -3 mm 2 s -1 for free water and d restricted = 1 × 10 -3 mm 2 s -1 for water diffusion restricted in cell bodies.

In addition to those isotropic water compartments, a number of anisotropic compartments (three in our experiments but again this number may be changed depending on the application) referred as DDI compartments, are added to model how water molecules diffuse in axonal bundles with various orientations. Diffusing water molecules in a particular anisotropic compartment are assumed to undergo a random displacement that is the independent sum of a von Mises-Fisher (VMF) random variable on the sphere S 2 of radius r and a Gaussian random variable in R 3 [START_REF] Stamm | A new multi-fiber model for low angular resolution diffusion MRI[END_REF]. The VMF is an orientation distribution with the following PDF on the unit sphere [START_REF] Banerjee | Clustering on the unit hypersphere using von Mises-Fisher distributions[END_REF]:

f (x|µ, κ) = κ 4π sinh(κ) exp(κµ T x) ( 2
)
where sinh is the hyperbolic sine, x ∈ S 2 , µ ∈ S 2 is the principal axis of diffusion, κ ∈ R + is an orientation dispersion index around µ (see Fig 1). This definition of a VMF distribution is extended to a sphere of radius r: Hence, the resulting PDF describing this random displacement is given by the 3D convolution of the VMF distribution with the Gaussian distribution:

f (x|µ, κ, r) = κ 4π sinh(κ)r 3 2 exp κµ x r (3)
p 0 (x|µ, κ, d, ν) = f (x|µ, κ, r) * exp -1 2 x T Σ -1 x (2π) 3 2 |Σ| 1 2 (4)
where r is the radius of the VMF sphere given by r = √ νd and Σ the covariance matrix of the Gaussian part defined as:

Σ = (1 -ν)d κ + 1 [I 3 + κµµ T ] ( 5
)
where d is the diffusivity along vectors µ and ν the non Gaussian part of the compartment. The VMF distribution has a dominant direction along µ. However an anisotropic compartment represents a bundle of axons that are assumed to be in a unique orientation with no privileged forward or backward direction within a voxel. Therefore the water diffusion is equal in both directions ±µ. To consider orientation instead of direction, the PDF is then symmetrized giving the final PDF of the DDI compartment [START_REF] Stamm | A new multi-fiber model for low angular resolution diffusion MRI[END_REF]:

p(x|µ, κ, d, ν) = p(x|µ, κ, d, ν) + p(x| -µ, κ, d, ν) 2 (6)

III. AVERAGING AND INTERPOLATION AS A SIMPLIFICATION PROBLEM

A. Global interpolation scheme

Interpolation or averaging of MCMs from different images may seem to be different problems. However they can all be treated as the computation of a MCM best matching a linear combination of MCMs. Let us consider two examples, one of each. If we consider the construction of an atlas, i.e. an MCM image summarizing a collection of K MCM images registered in a common space, the construction of the mean image requires to average at each voxel several MCMs (each from an image with a given weight, often 1/K but which may vary in the case of spatio-temporal atlases) into only one. The second example can be the registration of an MCM image onto another or the resampling of it by a transformation. In this case, we need to define an interpolation scheme for MCMs. Trilinear interpolation is the most classical approach in the literature for this. However this method, as well as other interpolation methods, ends up at each given voxel with a set of several MCMs (K = 8 for trilinear interpolation), each with a weight, that we need to summarize into a single MCM. The only difference with the averaging example is that the weights vary depending on the spatial location in the image. Since voxels can be treated separately in both tasks, interpolation and averaging for MCMs are in fact solved if we are able to define a scheme for the computation of a single MCM from multiple MCMs each with a weight.

We therefore consider the general problem of obtaining, from a linear combination of K MCMs (each weighted by a weight α k , with K k=1 α k = 1), an MCM that best summarizes the information contained in the combination. A simple solution, since MCMs are already a linear combination of individual compartments, would be to construct an MCM containing all compartments from the K models, each weighted according to the α k values and internal weights.

However, this is not satisfactory as this would result in an intractable model containing a large number of compartments, each not meaning much about the underlying microstructure of the brain. In this section, we therefore tackle the interpolation and averaging underlying problems as a simplification problem with the following aim: best summarizing the information in a linear combination of MCMs into a single MCM with few compartments (ideally the same number as the input MCMs but that parameter is left to the user).

We treat this problem by following the scheme illustrated in Fig. 2. This scheme is divided into two parts (treated in the following sub-sections and with a proposed implementation in Section IV-B for the studied models in the experiments): isotropic compartments simplification and anisotropic compartments simplification. The reason for this separation is the following. While anisotropic compartments share the same interpretation (a fiber bundle inside the voxel with given microstructural properties), isotropic compartments each describe a cell or tissue type (e.g. grey matter, glial cells, CSF...). We have therefore chosen not to mix their simplification with anisotropic compartments to keep the interpretability of the resulting model.

B. Isotropic compartments simplification

In addition and for the same reason as above, isotropic compartments are treated separately as they each model a specific tissue type and should not be mixed. For each of the M output isotropic compartments, we want to summarize the information coming from the K input MCMs into a single compartment of the same type. This requires 1-to compute the output weight ŵiso,l of each isotropic compartment in the final MCM, 2-to simplify a weighted sum of isotropic compartments into a single one.

The first part of the problem, weight computation, is relatively straightforward. It amounts, for the isotropic compartment of index l, to multiplying the input weights w iso,l of each input MCM by the relative weights of the K input MCMs:

ŵiso,l = K k=1 α k w iso,l,k (7)
where w iso,l,k denotes the weight w iso,l in Eq. ( 1) for the k-th MCM. With this done, the second part is model specific but only requires to know a metric space or technique to simplify M isotropic compartments into a single one. This second part is described in Section IV-A together with its simplification for the models considered in our experiments.

C. Anisotropic compartments simplification

Anisotropic compartments are more complex to deal with since there are several in each input MCM (from zero to usually three to model realistically crossing fibers in the brain), and their number may vary due to e.g. model selection in the MCM estimation process [START_REF] Stamm | Fast identification of optimal fascicle configurations from standard clinical diffusion MRI using Akaike information criterion[END_REF], [START_REF] Behrens | Characterization and propagation of uncertainty in diffusion-weighted MR imaging[END_REF]. For these compartments, our framework takes as an input a number N of output anisotropic compartments in the simplified MCM. This number, left to the user can e.g. be three by default to ensure we keep complex fiber crossings such as the ones in the centrum semiovale [START_REF] Budde | Quantification of anisotropy and fiber orientation in human brain histological sections[END_REF], [START_REF] Wedeen | The geometric structure of the brain fiber pathways[END_REF]. The problem here is therefore a bit different than for isotropic compartments: we start from a set of K MCMs, each with N (k) anisotropic compartments and wish to simplify them into N compartments.

Compartments clustering. This problem is a clustering problem. Reducing a number of objects or variables into a smaller group (often denoted as clusters) is a common issue to many different domains: machine learning, data compression or image segmentation for example [START_REF] Haralick | Image segmentation techniques[END_REF]. This is a difficult and open problem with a large literature dedicated to it [START_REF] Jain | Data clustering: a review[END_REF]. Here, we want to cluster a total of Q compartments (here

Q = k N (k)) into N compartments.
Among the most popular algorithms proposed, the k-means consists in minimizing the distance across groups between points and their cluster centroids [START_REF] Lloyd | Least squares quantization in PCM[END_REF] [START_REF] Likas | The global K-Means clustering algorithm[END_REF]. In this algorithm, each data point is affected to one and only one class. The algorithm thus alternates, from an initial position, between centroid computation and cluster affectation. The minimization of the distance in the k-means algorithm is a NP-hard problem [START_REF] Mahajan | The planar K-Means problem is NP-hard[END_REF], for which a large number of approximations in a polynomial time have been proposed [START_REF] Hartigan | Algorithm AS 136: A K-Means clustering algorithm[END_REF]- [START_REF] Park | A simple and fast algorithm for K-medoids clustering[END_REF]. However these methods are sensitive to a random initialization and the risk to converge to a local minimum is high [START_REF] Bradley | Refining initial points for K-Means clustering[END_REF]. Another clustering category considers a different way of affecting data points to classes by assigning a weight to each point which corresponds to a probability of membership to each cluster. The fuzzy C-means clustering belongs to this category [START_REF] Bezdek | FCM: The fuzzy C-Means clustering algorithm[END_REF], [START_REF] Pal | A possibilistic fuzzy C-Means clustering algorithm[END_REF].

"True" anisotropic MCM    s 11 • • • s 1n . . . . . . . . . s n1 • • • s nn    Similarity matrix [x 1 , ..., x n ] ∈ R n×k

Largest eigenvectors of the Laplacian

Simplification of the input PDF parameters according to the semimetric chosen for S Final weights computation: ŵl in Eq. ( 9)

β j,k,l
Weights assigned using fuzzy C-means Here, we use a modified version of normalized spectral clustering algorithm [START_REF] Ng | On spectral clustering: Analysis and an algorithm[END_REF], illustrated in Fig. 3. From a similarity matrix S between inputs, spectral vectors of dimension N are computed and used to initialize a clustering algorithm. Here, we choose a fuzzy C-means, instead of the regular k-means algorithm, to keep a continuity when input data weights vary. Several choices are possible to define spectral vectors : we use the largest eigenvectors of the symmetric normalized Laplacian of S [START_REF] Ng | On spectral clustering: Analysis and an algorithm[END_REF]. This method only needs a similarity matrix between the inputs (here MCM compartments), is robust, well-reviewed with a large number of algorithms proposed [START_REF] Luxburg | A tutorial on spectral clustering[END_REF], [START_REF] Spielmat | Spectral partitioning works: Planar graphs and finite element meshes[END_REF]. Considering a specific distance (or semi-metric depending on the compartment type) between two anisotropic compartments PDFs {d(p i , p j )} i,j=1,...,Q , the similarity matrix S is computed using the following:

S i,j = exp -d 2 (p i , p j )) 2σ 2 i,j=1,...,Q (8)
where d has at least the same properties than classic distance without the triangle inequality. This defines a weak metric in spaces not equipped with a natural distance. σ is a normalization coefficient to avoid computational errors. From this similarity matrix, the fuzzy C-means algorithm is able provide us with membership weights, describing how much j-th compartment from the k-th input MCM fits into a l-th output class. We denote these membership weights by β j,k,l .

Final weights computation With these weights defined, and similarly to the previous section, the simplification problem is solved in two parts: 1-output weights ŵl computation (for l = 1 to N ) and output models parameters computation for each output anisotropic compartment PDF pl . The weights are computed as follows:

ŵl = K k=1 α k N (k) j=1 w j,k β j,k,l (9)
where w j,k is the weight w j in Eq. ( 1) for the k-th MCM.

With the additional constraints that ∀k = 1..K, ∀j = 1..N (k),, l β j,k,l = 1 (common to all clustering algorithms), that the initial weights of each input MCM sum up to 1 ( i w iso,i,k + j w j,k = 1), and that k α k = 1, we can verify that the output weights (isotropic and anisotropic) sum up to 1. Again, the second part of the problem is the estimation of the N output anisotropic compartments i.e. the pl PDFs parameters. This part relies on the definition of a semi-metric between compartments PDFs and the definition of an associated weighted mean. A solution to this problem is derived for the two studied models in the next sections: the tensor compartment in Section IV-B.1 and the DDI compartment in Section IV-B.2.

D. Summary of required bricks for simplification

This global simplification process may seem complex to tackle. In fact, it requires only few pieces to be defined to be implemented for virtually any model and is then very versatile. This is thanks to the spectral clustering approach chosen above which, contrarily to Taquet et al. approach [START_REF] Taquet | Registration and analysis of white matter group differences with a multi-fiber model[END_REF], does not assume the PDF of the anisotropic compartments or isotropic compartments to belong to the exponential family. We end this global framework presentation by recalling the few bricks needed to implement this simplification approach:

• A distance or semi-metric between isotropic compartments of each type (usually simple as most of these compartments are in fact Gaussian PDFs) • A distance or semi-metric between two anisotropic comd(p i , p j ) • A way to compute resulting PDF parameters from a set of weighted input PDFs for both isotropic and anisotropic

IV. COMPARTMENTS SIMPLIFICATION IMPLEMENTATION

A. Isotropic compartments

Isotropic compartments may take on different forms [START_REF] Panagiotaki | Compartment models of the diffusion MR signal in brain white matter: a taxonomy and comparison[END_REF]. With the exception of a few (e.g. the Stanisz compartment [START_REF] Stanisz | An analytical model of restricted diffusion in bovine optic nerve[END_REF]), most of these compartments follow a zero-mean isotropic 3D Gaussian distribution:

p iso (x) = 1 (2πd iso ) 3 2 exp - x x 2d iso ( 10
)
where d iso is the diffusivity. In our method implementation, we will consider only this kind of Gaussian isotropic compartments.

We therefore want, for each tissue type, to simplify a set of isotropic zero-mean Gaussian PDFs into a single one. This falls back in the framework of symmetric positive definite matrices processing, for which a number of Riemannian or Lie group approaches have been proposed. We have chosen in this paper to rely on the log-Euclidean framework [START_REF] Arsigny | Log-Euclidean metrics for fast and simple calculus on diffusion tensors[END_REF] for its fast computation capabilities. The computations are even simpler for isotropic Gaussian PDFs and amount to a geometric mean of the input diffusivities:

diso,l = exp 1 ŵiso,l K i=1 α k w iso,l,k log(d iso,l,k ) (11)
with d iso,l,k is the isotropic diffusivity of the l-th isotropic compartment of the k-th MCM. Additionally, one may note that in the case of fixed, equal, diffusivities for a given tissue type as it may happen in practice especially for models estimated on clinical datasets, the output diffusivity will be the same as the input diffusivities and only the weight of the given compartment will need to be computed.

B. Anisotropic compartments

Each MCM has one or several anisotropic compartments that represent a constrained water diffusion along particular tissues. Within an axon, the water is indeed trapped and then diffuses along the axon direction. A bundle of axons oriented in the same direction creates an anisotropic water diffusion at the voxel level. This water diffusion is represented as an anisotropic compartment following different models depending on the MCM considered. We present here the distances or semi-metrics, together with the computation of the output compartment PDFs parameters for the tensor compartment (that gives rise to the multi-tensor model -MTM) and the diffusion directions imaging compartment (that gives rise to the DDI model).

0278-0062 (c) 2020 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information. For each cluster k, we wish to summarize the set of PDFs p j,k (PDFs of the j-th compartment in the k-th input MCM) with their associated weights from Eq. ( 9) into a single compartment of the same type. To simplify notations in the following section, we now just consider n anisotropic compartments with their corresponding weights w i that we want to simplify into a single one. To perform spectral clustering, we also need to be able to compute a distance or semi-metric between two anisotropic compartments. We do so in the following sub-sections for the two compartment types.

1) Tensor compartments: The MTM is a straightforward extension of the classic DTI. In addition to isotropic compartments, the anisotropic compartments are modeled as tensors, i.e, 3D Gaussian distributions with zero-mean, parameterized by a symmetric positive-definite matrix. Let {T i } i=1..n be a set of tensors with their corresponding weights {w i } i=1..n . Estimating the output PDF parameters amounts here to compute a tensor T from this set of tensors and weights. All {T i } i=1..n belong to S + 3 (R), the space of positive-definite matrices. On this space, we can use again the log-Euclidean framework [START_REF] Arsigny | Log-Euclidean metrics for fast and simple calculus on diffusion tensors[END_REF]. The simplified tensor is recovered as:

T = exp n i=1 w i log(T i ) (12)
From this, we can easily and efficiently simplify any number of tensors into a unique one. A distance between two tensors is also directly defined in the log-Euclidean space as:

d(T 1 , T 2 ) = || log(T 1 ) -log(T 2 )|| F (13)
where ||.|| F is the Frobenius norm. The similarity matrix used for spectral clustering is thus derived from Eq. 8.

2) Diffusion direction imaging compartments: Again, to simplify notations, we consider n DDI compartments with their corresponding weights w i that we want to simplify into one compartment. Let µ i , κ i , ν i , d i be the parameters of the input compartments (refer to Section II-C for more details) and µ, κ, ν, d be the parameters of the final compartment.

This distribution is too complex to allow for a natural analytic metric. Thus, we propose 4 different methods to perform the simplification based on the different parts of the DDI compartment PDF: (1) simplest, a naive method using the Euclidean metric; (2) tensor, a metric based on the compartments main orientations; (3) log VMF, a metric based on the VMF distribution, (4) covariance analytic, a metric based on the Gaussian part of the DDI compartment. Simplest. Each µ i is a unit direction in S 2 . However, they do not represent a direction but an orientation. The simplest way to solve this problem (as two opposite directions) is to put all µ i in the top hemisphere and compute their Karcher mean on the sphere to obtain µ. The rest of the parameters are estimated as a weighted Euclidean mean.

Here, all parameters are estimated independently and we wish to compute a similarity measure to this semi-metric. A weighted sum of distances using all parameters is used and defined as:

d simple (p 1 , p 2 ) =| < µ 1 , µ 2 > | + δ|κ 1 -κ 2 |+ β|d 1 -d 2 | + γ|ν 1 -ν 2 | (14)
with δ, β, γ some terms to normalize the influence of the parameters in d simple . The similarity matrix S used for spectral clustering related to this approach is then defined using the semi-metric defined in Eq. ( 14) along Eq. ( 8).

This method should be considered only as the most naive approach to create a semi-metric and give a reference. In a general case, it does not constitute a good mathematical framework to tackle this problem.

Tensor. The simplest approach is only a partial solution, especially for directions close to the sphere equator, which might generate discontinuities. We now consider µ i as orientations instead of directions. To do so, µ i is represented as a cigar-shaped tensor T i defined as:

T i = µ i µ T i + εI 3 (15)
with ε = 10 -6 to have non degenerated tensors and I 3 the identity matrix. Then, T i ∈ S + 3 (R) and we can use the log-Euclidean framework on tensors. T is computed in a similar way to MTM (see Eq. 12). µ thus becomes the principal direction of T (i.e the eigenvector with the largest eigenvalue). The other parameters are obtained by the weighted Euclidean mean as for the simplest averaging. Also the semi-metric between compartments is defined as:

d tensor (p 1 , p 2 ) = || log(T 1 ) -log(T 2 )|| F + δ|κ 1 -κ 2 |+ β|d 1 -d 2 | + γ|ν 1 -ν 2 | (16)
where δ, β and γ are normalization coefficients.

Covariance analytic. Another approach uses information from covariance matrices Σ i of the DDI compartments. These Σ i matrices belong to S + 3 (R) and can be merged into Σ using the log-Euclidean framework:

Σ = exp n i=1 w i log(Σ i ) (17)
We then wish to extract all parameters from Σ. We start by approximating Σ by a cigar-shaped tensor to match the DDI compartment structure. To do this, we need to enforce two equal secondary eigenvalues λ ⊥ . In the log-Euclidean framework, this amounts to compute λ ⊥ as λ ⊥ = √ λ 2 λ 3 where λ 2 , λ 3 are the two lowest eigenvalues of Σ. We now have Σ the cigar-shaped tensor of C:

Σ = (1 -ν)d κ + 1 [I 3 + κµµ T ] (18)
A tensor in S + 3 (R) with 3 different eigenvalues has 6 specific parameters, however our tensor is cigar-shaped and then only has 4 specific parameters. We want to recover 2 parameters for the direction µ and 1 for κ, µ and d. For the missing 5-th 0278-0062 (c) 2020 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2020.3042765, IEEE Transactions on Medical Imaging parameter, we process by identification using the relationship r 2 = νd. To exploit this relationship, we define the weighted mean radius as the one whose sphere surface is the weighted mean of the input sphere surfaces. This corresponds to a weighted Euclidean mean of the individual r 2 i :

r = n i=1 w i r 2 i 1 2 (19)
This therefore gives us a direct relation between ν and d leading to only 4 parameters to estimate (µ, κ and ν), d being computed as d = r 2 /ν. Then, we can estimate all the parameters by the resolution of the eigensystem of Σ. We obtain µ as the principal eigenvector of Σ and the following estimation for ν and κ :

     ν = r 2 λ + r 2 κ = λ λ -λ ⊥ (20)
In the covariance analytic method, all the parameters except r are estimated from the eigen analysis of Σ, the weighted mean of the covariance matrices. Thus, the natural choice to define a distance in this case is to compute directly the log-Euclidean distance:

d covariance analytic (p 1 , p 2 ) = || log(Σ 1 ) -log(Σ 2 )|| F (21)
log von Mises Fisher. We now explore the option to use the VMF to compute µ and κ and recover only ν from Σ. We want to consider a VMF distribution as a point in a Riemannian manifold following geometric information methods. To define a metric on PDF spaces, a theoretical framework exists using partial derivatives of the PDF though Christoffel symbols [START_REF] Amari | Methods of information geometry[END_REF]. We tried to apply such a framework to VMF distributions, unfortunately, the corresponding partial derivative equations are not analytically solvable. Therefore, to consider a VMF as a point of a Riemannian manifold, the two parameters need to be treated separately. The space S 2 × R + , where the couple of parameters (µ, κ) evolves, is a Riemannian manifold as a product of two Riemannian manifolds. Therefore the averaging of the two parameters µ and κ can be treated independently. To compute the weighted mean of several points, a geodesic on these manifolds is defined (refer to [START_REF] Mcgraw | von Mises-Fisher mixture model of the diffusion ODF[END_REF] for details). Following this geodesic, the orientation averaging is similar to tensor averaging. The interpolation of κ is done recursively by projection as in McGraw et al. [START_REF] Mcgraw | von Mises-Fisher mixture model of the diffusion ODF[END_REF]. Letting κ = κ 1 , we repeat until convergence (i.e until l κ < ε):

     l κ = n i=1 w i log κ i κ κ = κ exp(l κ ) (22)
This iterative method converges to the Fréchet mean in the tangent space. Similarly to the covariance analytic method, we then use the relation r 2 = νd. Then knowing all parameters except ν, we obtain it from Σ as in the covariance analytic method. As 2 equations are available to estimate one parameter, we compute the least square solution to ensure more stability:

ν = r 2 [2r 2 + λ + λ ⊥ (1 + κ)] 2(r 2 + λ)[r 2 + λ ⊥ (1 + κ)] (23)
The log VMF metric in the Riemannian manifold S 2 × R + belongs to the log-space for both direction and concentration parameters. Therefore we compute the distance of tensors and the distance of κ with their corresponding Riemannian metric which give the following semimetric between compartments:

d logVMF (p 1 , p 2 ) = || log(T 1 ) -log(T 2 )|| F + δ| log(κ 1 ) -log(κ 2 )| + β||r 1 -r 2 || 2 2 (24)
where δ, β and γ are normalization coefficients.

V. EXPERIMENTS AND RESULTS

A. Evaluation of DDI compartment simplification on simulated data

The MCM simplification holds two different parts: the global scheme which is generic and applicable to any MCM and the compartments simplification part that is related to the anisotropic compartment considered. The evaluation of the global interpolation for two MCM, the MTM and the DDI, is done in Section V-B. In this section, we present the evaluation of the compartments simplification on the DDI compartment according to the different metrics proposed.

To do so, we evaluate the simplification 4 weighted DDI compartments into a single one on a 11 × 11 grid, illustrated in Fig. 4. The 4 corners of this grid correspond to 4 different reference DDI compartments. The rest of the grid is composed with the 4 DDI compartments of the corners weighted by their distance to the current spatial position (see reference in Fig 4). The orientation of each DDI compartment in the corners is randomly chosen on the sphere and the rest of the parameters are picked within the following uniform distributions: [0, 1] for ν, [0, 20] for κ, [5.10 -4 , 5.10 -3 ] mm.s -2 for d. From this reference grid including 4 compartments in each pixel, we compute for each method proposed a simplified grid with only one compartment: simplest, tensor, logVMF and covariance analytic. To perform a quantitative evaluation, in each pixel, a DWI having 180 directions on 3 shells (b = 1000, 2000, 3000 s.mm -2 ) is simulated for the reference and each evaluated simplified model following Eq. ( 6) in [START_REF] Stamm | A new multi-fiber model for low angular resolution diffusion MRI[END_REF]. Finally, the Euclidean error map is computed between the reference and each metric.

On the grids presented in Fig. 4, compartment orientations obtained by the simplest method are clearly not in accordance with the ground truth leading to large Euclidean errors. The orientations obtained by the tensor and logVMF methods seem more accurate, although they are still leading to large residual errors. In particular for these methods, regions with orthogonal directions keep large κ values which is not realistic. Crossing fibers should indeed induce more dispersion, which is more equivalent to a single compartment with a low κ value and a low d value. The covariance analytic method reflects this expected behavior resulting in a smaller error than all other methods. We further simulated 500 random grids and computed their simplification. Again, the Euclidean difference between the simulated DWI from these maps was computed, averaged and normalized so that the simplest method gets a reference result of 100. We obtain the following normalized values for the 4 methods: simplest 100, tensors 31.6, logVMF 28.0, covariance analytic 11.1. These results confirm that the covariance analytic method obtains the lowest error for the compartment simplification part.

B. Experiments on real data for MCM interpolation

1) MTM: now test the entire MCM interpolation pipeline including anisotropic and isotropic simplification. A MTM is estimated from a subject of the HCP data [START_REF] Van Essen | The Human Connectome Project: a data acquisition perspective[END_REF], which is a DWI with 145 × 174 × 145 voxels with a 1.25 × 1.25 × 1.25 mm 3 resolution and 270 gradient directions over 3 b-values (1000, 2000, and 3000 s.mm -2 ). The estimated MTM includes 3 anisotropic compartments (tensors) and two isotropic compartments (one free water with a diffusivity d free = 3×10 -3 mm 2 s -1 and one restricted water with a diffusivity d restricted = 1 × 10 -3 mm 2 s -1 ).

To test the global interpolation scheme, a rotation of angle 120 degrees is picked around a random axis and then applied 3 consecutive times to the original MTM image. We then compare the final MTM obtained to the original one. To visualize the MTM, we compute the MTM PDF values on several points of a sphere and deform it using these values. As a consequence, the sphere will be elongated along the most probable diffusion directions and contracted elsewhere. To perform a quantitative evaluation, the DWI corresponding to the original and rotated MCM were simulated using 180 gradient directions with 3 shells (b = 1000, 2000, 3000 s.mm -2 ). The absolute difference between the two simulated DWI were then averaged along the 180 gradients. The visual representation of the original and interpolated MTM is presented in Fig. 5. At the brain level, the MTM seem very similar to the original one although smoother. In the zoomed area, despite crossing in the original MTM, the rotated MTM stay close to the original models even in the crossing region. As expected, the image is smoothed by the interpolation but all main orientations are recovered. Visually after 3 consecutive interpolations the result seems very good compared to the original MTM. Quantitative results for this experiment are presented in Fig. 6. The DWI absolute difference is more important within the ventricles, which have an isotropic diffusion along with a high noise level. Yet, within the rest of the brain, the absolute difference is mostly under 10% which seems reasonable.

2) DDI: We then tested the entire MCM interpolation pipeline for the DDI. To perform the validation of the different DDI compartment interpolation methods on real data, we tested methods on a set of 46 real DDI estimated from DWI with 128×128×55 voxels with a 2×2×2 mm 3 resolution, 30 gradient directions with one b-value = 1000 s.mm -2 . Input DDI have been estimated with 3 DDI compartments and one free water compartment [START_REF] Stamm | Comprehensive maximum likelihood estimation of diffusion compartment models towards reliable mapping of brain microstructure[END_REF].

For each input DDI, we compute a rotation of angle 120 degrees and then apply it 3 consecutive times for each 4 The logVMF suffers from the multiple interpolations and obtains worse results than the simplest method, which stays far from the two best methods. The covariance analytic performs significantly better than the tensor showing an improvement, moderate but very robust, for several successive interpolations. Following the same process than for MTM, we visualize the original DDI and the one rotated, for one of the 46 subjects of our dataset, by the two methods that obtained the best results, tensor and covariance analytic(see Fig. 9). On the first line, the 3 images have the same aspect though the rotated DDI seems smoother than the original one, however more interesting observations come with a zoom. On the second line, the covariance analytic method compartments are smaller than the original ones probably due to an under estimation of the diffusivity. On the same line, the compartments size of the tensor method seem slightly bigger than the original DDI. The third line represents all the compartments at the same size to focus on the orientation. We can see that orientations of the compartments obtained with covariance analytic methods are very similar to the original ones while the compartments from the tensor method orientations seem different. These observations show that even if all these methods have imperfections, the covariance analytic method is able to robustly recover the orientation after 3 successive rotations.

C. DDI Atlas Construction

One of the goals of MCM registration is the production of an average atlas of the white matter microstructure. We computed an atlas from 46 DDI following Guimond et al. atlas construction method [START_REF] Guimond | Average brain models: A convergence study[END_REF]. This atlas construction was performed using non linear DTI registration as proposed by Suarez et al. [START_REF] Suarez | Automated delineation of white matter fiber tracts with a multiple region-of-interest approach[END_REF]. Then, the obtained transformations were applied to the DDIs. We interpolated the DDIs using our clustering method with the covariance analytic approach. In addition, when applying a transformation to oriented models, it is necessary to apply the local linear part of the transformation to the interpolated models. We used a technique similar to finite-strain reorientation for tensors [START_REF] Ruiz-Alzola | Nonrigid registration of 3D tensor medical data[END_REF] by applying the local rotation to the µ i directions of each anisotropic compartment of the interpolated DDI. We present the visual result of the atlas and a zoomed area in Fig. 10. This atlas provides a clear distinction of crossing fibers and will be of great interest in future studies for example of white matter microstructure destruction in diseases.

VI. DISCUSSION

In this study, we presented a general framework to interpolate MCM and we applied it for two anisotropic compartments: the tensor compartment and the DDI compartment. The averaging of several tensors can be done in a multiple ways preserving different properties of the tensors [START_REF] Wang | Tracking on the product manifold of shape and orientation for tractography from diffusion MRI[END_REF], [START_REF] Kindlmann | Geodesic-loxodromes for diffusion tensor interpolation and difference measurement[END_REF]. We chose the log-Euclidean framework as a fast and efficient method and more importantly because it offers a natural metric. On the contrary, no natural metric is available yet for the DDI compartment and it is more challenging to find a robust one as demonstrated in this paper. Thus, a potential direction for future work could be to investigate this question for other anisotropic compartments (NODDI, CHARMED, ...).

One other important point in our approach is the choice of the semi-metric or metric that is applied. This choice is crucial as it may change the parameters of the simplified models. It therefore has to be adapted to the application for which the interpolation or averaging is made. We have chosen in our experiments and method to use metrics that preserves as much as possible the PDF of the resulting simplified compartments but other options could be explored such as more microstructure oriented metrics that would preserve other properties of the compartments.

Another acute problem is the number of clusters of the model. In our experiments, the number of anisotropic compartments is the same for the entire image. A priori, some parts of the brain do not need 3 compartments in addition to a free water compartment. Tools based on AIC exist to compute a specific number of compartments for each voxel [START_REF] Sakamoto | Akaike information criterion statistics[END_REF]. The number of compartments after interpolation in our algorithm is however fixed for the entire image. Simple solutions towards handling better images with varying number of anisotropic compartments include taking the mean or maximum of the number of clusters of the input voxels do not support successive transformations. We can imagine smart algorithms to estimate automatically the optimal number of clusters after interpolation. Such methods already exist for spectral clustering based on rotation of the spectral vector [START_REF] Zelnik-Manor | Self-tuning spectral clustering[END_REF] or the selection of eigenvalues [START_REF] Sanguinetti | Automatic determination of the number of clusters using spectral algorithms[END_REF]. Unfortunately they cannot consider one cluster as a solution and are thus not directly applicable to our problem.

Finally, the purpose of this method is to provide efficient tools to directly register a set of MCMs on a common reference. A clinical study often includes data from a group of healthy subjects, and a group of patients affected by the same disease. From such a database, we are looking for biomarkers that highlight anatomical abnormalities. The registration of such abnormalities, as lesions, is an important topic and difficult one. There is indeed no guarantee that a lesion after the registration process can be detectable or keep the same properties. This complex problem is common to all the registration methods, including ours. Therefore, it could be interesting, as a future work, to compare how the preservation of lesion biomarkers depends on the chosen registration method.

VII. CONCLUSION

As MCM becomes increasingly used in research studies and are on the line to be evaluated on clinical one, tools specifically designed for such models are required. In this paper, we proposed a method to interpolate and average MCM images as a simplification problem. Isotropic and anisotropic compartments averaging are treated separately, this latter part using spectral clustering adapted for this purpose. This simplification process, very generic, can be applied to any MCM provided that we are able to compute semi-metrics between anisotropic compartments. The code for MCM image averaging and interpolation is available open-source as part of the Anima 1 software package.

We tested the entire pipeline applying 3 consecutive rotations to an original MTM estimated from HCP data. For each applied rotation, at each voxel, 24 anisotropic compartments are simplified into 3 ones using spectral clustering. Yet, after 3 rotations the interpolated MTM is very similar compared to the original MTM showing its ability to preserve crossing fibers. In addition to this, quantitative experiments were also 1 Anima: Open source software for medical image processing from the Empenn team. https://anima.irisa.fr -RRID:SCR 017017 performed on DDI comparing 4 different metrics for both synthetic and real data. The covariance analytic method exhibits significantly better performance than other ones for both experiments. As we saw in the visual analysis there are still some questions to be analyzed further on the attenuation of the diffusivity with the covariance analytic method, nevertheless it is able to robustly recover the orientation after several transformations.

0278-0062 (c) 2020 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission. See http://www.ieee.org/publications_standards/publications/rights/index.html for more information.

Fig. 1 :

 1 Fig. 1: Points sampled from 3 VMF distributions on the sphere (blue: κ = 200, green: κ = 20, red: κ = 2). The mean directions µ are shown with arrows.

Fig. 2 :

 2 Fig. 2: Global scheme of the simplification of 4 MCMs into a single output MCM with a two steps process: isotropic simplification and anisotropic simplification.

Fig. 3 :

 3 Fig. 3: Illustration of spectral clustering for anisotropic compartments simplification.

  This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2020.3042765, IEEE Transactions on Medical Imaging IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. XX, NO. XX, XXXX 2020

Fig. 4 :

 4 Fig. 4: Illustration of DDI compartment simplification using the 4 methods superimposed on their local error maps with the corresponding scalar bar at the bottom right. The original grid with 4 superposed compartments at each pixel is the reference image.

Fig. 5 :

 5 Fig. 5: Visualization of an original MTM and its corresponding interpolation after 3 rotations.

Fig. 6 :

 6 Fig. 6: First column: original MTM; second column: interpolated MTM after 3 rotations of 120 degrees; third column: absolute difference between the DWI simulated from the original and rotated MTM.

Fig.

  Fig. Illustration of DDI interpolation with the 4 proposed semi-metrics. The evaluation is done after 3 consecutive rotations of 120 degrees. Each method is represented by two contiguous images, the interpolated DDI and the DWI Euclidean distances. The original DDI is located at the center.

Fig. 8 :

 8 Fig. 8: Error between the DWI corresponding images of the original DDI and the one estimated after 3 rotations following the 4 different methods.

Fig. 9 :

 9 Fig. 9: Illustration of two DDI compartment interpolation methods, covariance analytic and tensor, compared to the original DDI. The interpolation is made applying 3 consecutive rotation of 120 degrees. The first line represent the original DDI and the two interpolated by the covariance analytic method and the tensor method. The second line represents the same images zoomed on the yellow rectangle. The third line represents the same zoom with all the compartments normalized to the same size to focus on the orientation of the DDI compartments.

Fig. 10 :

 10 Fig. 10: Example of a DDI atlas superimposed on the average B0 image: Axial view on the first line and coronal view. This atlas is constructed from our database of 46 real DDI images.

  This article has been accepted for publication in a future issue of this journal, but has not been fully edited. Content may change prior to final publication. Citation information: DOI 10.1109/TMI.2020.3042765, IEEE Transactions on Medical Imaging IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. XX, NO. XX, XXXX 2020

Authorized licensed use limited to: INRIA. Downloaded on January 14,2021 at 15:24:36 UTC from IEEE Xplore. Restrictions apply.

IEEE TRANSACTIONS ON MEDICAL IMAGING, VOL. XX, NO. XX, XXXX 2020

This manuscript is submitted for review on Renaud H édouin PhD is partially funded by the Brittany Council. MRI data acquisition was supported by the Neurinfo MRI research facility from the University of Rennes I. Neurinfo is granted by the the European Union (FEDER), the French State, the Brittany Council, Rennes Metropole,