Xavier Pennec

Pierre Fillard

Nicholas Ayache

A Riemannian Framework for Tensor Computing

Keywords: Tensors, Diffusion Tensor MRI, Regularization, Interpolation, Extrapolation, PDE, Riemannian Manifold, Affine-invariant Metric

Tensors are nowadays a common source of geometric information. In this paper, we propose to endow the tensor space with an affine-invariant Riemannian metric. We demonstrate that it leads to strong theoretical properties: the cone of positive definite symmetric matrices is replaced by a regular and complete manifold without boundaries (null eigenvalues are at the infinity), the geodesic between two tensors and the mean of a set of tensors are uniquely defined, etc. We have previously shown that the Riemannian metric provides a powerful framework for generalizing statistics to manifolds. In this paper, we show that it is also possible to generalize to tensor fields many important geometric data processing algorithms such as interpolation, filtering, diffusion and restoration of missing data. For instance, most interpolation and Gaussian filtering schemes can be tackled efficiently through a weighted mean computation. Linear and anisotropic diffusion schemes can be adapted to our Riemannian framework, through partial differential evolution equations, provided that the metric of the tensor space is taken into account. For that purpose, we provide intrinsic numerical schemes to compute the gradient and Laplace-Beltrami operators. Finally, to enforce the fidelity to the data (either sparsely distributed tensors or complete tensors fields) we propose least-squares criteria based on our invariant Riemannian distance which are particularly simple and efficient to solve.

Introduction

Positive definite symmetric matrices (so-called tensors in this article) are often encountered in image processing, for instance as covariance matrices for characterizing statistics on deformations, or as an encoding of the principal diffusion directions in Diffusion Tensor Imaging (DTI). The measurements of these tensors is often noisy in real applications and we would like to perform estimation, smoothing and interpolation of fields of this type of features. The main problem is that the tensor space is a manifold that is not a vector space with the usual additive structure. As symmetric positive definite matrices constitute a convex half-cone in the vector space of matrices, many usual operations (like the mean) are stable in this space. However, problems arise when estimating tensors from data (in standard DTI, the estimated symmetric matrix could have negative eigenvalues), or when smoothing fields of tensors: the numerical schemes used to solve the Partial Differential Equation (PDE) may sometimes lead to negative eigenvalues if the time step is not small enough. Even when a SVD is performed to smooth independently the rotation (eigenvectors basis trihedron) and eigenvalues, there is a continuity problem around equal eigenvalues.

In previous works [Pennec, 1996, Pennec andAyache, 1998], we used invariance requirements to develop some basic probability tools on transformation groups and homogeneous manifolds. This statistical framework was then reorganized and extended in [Pennec, 1999, Pennec, 2004] for general Riemannian manifolds, invariance properties leading in some case to a natural choice for the metric. In this paper, we show how this theory can be applied to tensors, leading to a new intrinsic computing framework for these geometric features with many important theoretical properties as well as practical computing properties.

In the remaining of this section, we quickly investigate some connected works on tensors. Then, we summarize in Section 2 the main ideas of the statistical framework we developed on Riemannian manifolds. The aim is to exemplify the fact that choosing a Riemannian metric "automatically" determines a powerful framework to work on the manifold through the introduction of a few tools from differential geometry. In order to use this Riemannian framework on our tensor manifold, we propose in Section 3 an affine-invariant Riemannian metric on tensors. We demonstrate that it leads to very strong theoretical properties, as well as some important practical algorithms such as an intrinsic geodesic gradient descent. Section 4 focuses on the application of this framework to an important geometric data processing problem: interpolation of tensor values. We show that this problem can be tackled efficiently through a weighted mean optimization. However, if weights are easy to define for regularly sampled tensors (e.g. for linear or tri-linear interpolation), the problem proved to be more difficult for irregularly sampled values.

With Section 5, we turn to tensors field computing, and more particularly filtering. If the Gaussian filtering may still be defined through weighted means, the partial differential equation (PDE) approach is slightly more complex. In particular, the metric of the tensor space has to be taken into account when computing the magnitude of the spatial gradient of the tensor field. Thanks to our Riemannian framework, we propose efficient numerical schemes to compute the gradient, its amplitude, and the Laplace-Beltrami operator used in linear diffusion. We also propose an adjustment of this manifold Laplacian that realizes an anisotropic filtering. Finally, Section 6 focuses on simple statistical approaches to regularize and restore missing values in tensor fields. Here, the use of the Riemannian distance inherited from the chosen metric is fundamental to define least-squares data attachment criteria for dense and sparsely distributed tensor fields that lead to simple implementation schemes in our intrinsic computing framework.

Related work

Quite an impressive literature has now been issued on the estimation and regularization of tensor fields, especially in the context of Diffusion Tensor Imaging (DTI) [START_REF] Basser | MR diffusion tensor spectroscopy and imaging[END_REF], Le Bihan et al., 2001, Westin et al., 2002]. Most of the works dealing with the geometric nature of the tensors has been performed for the discontinuity-preserving regularization of the tensor fields using Partial Differential Equations (PDEs). For instance, [START_REF] Coulon | Diffusion tensor magnetic resonance image regularization[END_REF] anisotropically restores the principal direction of the tensor, and uses this regularized directions map as an input for the anisotropic regularization of the eigenvalues. A quite similar idea is adopted in [Tschumperlé, 2002], where a spectral decomposition W (x) = U (x) D(x) U (x) T of the tensor field is performed at each points to independently regularize the eigenvalues and eigenvectors (orientations). This approach requires an additional reorientation step of the rotation matrices due to the non-uniqueness of the decomposition (each eigenvector is defined up its sign and there may be joint permutations of the eigenvectors and eigenvalues) in order to avoid the creation of artificial discontinuities. Another problem arises when two or more eigenvalues become equal: a whole subspace of unit eigenvectors is possible, and even a re-orientation becomes difficult. An intrinsic integration scheme for PDEs that uses the exponential map has been added in [START_REF] Chefd | Constrained flows of matrix-valued functions: Application to diffusion tensor regularization[END_REF], and allows to perform PDEs evolution on the considered manifold without re-projections. In essence, this is an infinitesimal version of the intrinsic gradient descent technique on manifolds we introduced in [Pennec, 1996, Pennec, 1999] for the computation of the mean.

The affine-invariant Riemannian metric we detail in Section 3.3 may be traced back to the work of [Nomizu, 1954] on affine invariant connections on homogeneous spaces. It is implicitly hidden under very general theorems on symmetric spaces in many differential geometry textbooks [START_REF] Kobayashi | Foundations of differential geometry[END_REF], Helgason, 1978, Gamkrelidze, 1991] and sometimes considered as a well known result as in [Bhatia, 2003]. In statistics, it has been introduced as the Fisher information metric [Skovgaard, 1984] to model the geometry of the multivariate normal family. The idea of the invariant metric came to the mind of the first author during the IPMI conference in 2001 [START_REF] Coulon | A regularization scheme for diffusion tensor magnetic resonance images[END_REF], Batchelor et al., 2001], as an application to diffusion tensor imaging (DTI) of the statistical methodology on Riemannian manifolds previously developed (and summarized in the next Section). However, this idea was not exploited until the end of 2003, when the visit of P. Thompson (UCLA, USA) raised the need to interpolate tensors that represent the variability from specific locations on sulci to the whole volume. The expertise of the second author on DTI [START_REF] Fillard | Quantitative analysis of white matter fiber properties along geodesic paths[END_REF] provided an ideal alternative application field. During the writing of this paper, we discovered that the invariant metric has been independently proposed by [START_REF] Förstner | A metric for covariance matrices[END_REF] to deal with covariance matrices, and very recently by [START_REF] Fletcher | Principal geodesic analysis on symmetric spaces: Statistics of diffusion tensors[END_REF] for the analysis of principal modes of sets of diffusion tensors. By looking for a suitable metric on the space of Gaussian distributions for the segmentation of diffusion tensor images, [START_REF] Lenglet | Statistics on multivariate normal distributions: A geometric approach and its application to diffusion tensor MRI[END_REF][START_REF] Lenglet | Toward segmentation of 3D probability density fields by surface evolution: Application to diffusion MRI[END_REF]] also end-up with the same metric. It is interesting to see that completely different approaches, relying on an affine-invariant requirement on the one hand, and relying on an information measure to evaluate the distance between distributions on the other hand, lead to the same metric on the tensor space. However, to our knowledge, this Riemannian metric has not been promoted as a complete computing framework, as we propose in this paper.

Statistics on geometric features

We summarize in this Section the theory of statistics on Riemannian manifolds developed in [Pennec, 1999, Pennec, 2004]. The aim is to exemplify the fact that choosing a Riemannian metric "automatically" determines a powerful framework to work on the manifold through the use of a few tools from differential geometry.

In the geometric framework, one can specify the structure of a manifold M by a Riemannian metric. This is a continuous collection of scalar products on the tangent space at each point of the manifold. Thus, if we consider a curve on the manifold, we can compute at each point its instantaneous speed vector and its norm, the instantaneous speed. To compute the length of the curve, we can proceed as usual by integrating this value along the curve. The distance between two points of a connected Riemannian manifold is the minimum length among the curves joining these points. The curves realizing this minimum for any two points of the manifold are called geodesics. The calculus of variations shows that geodesics are the solutions of a system of second order differential equations depending on the Riemannian metric. In the following, we assume that the manifold is geodesically complete, i.e. that the definition domain of all geodesics can be extended to R. This means that the manifold has no boundary nor any singular point that we can reach in a finite time. As an important consequence, the Hopf-Rinow-De Rham theorem states that there always exists at least one minimizing geodesic between any two points of the manifold (i.e. whose length is the distance between the two points).

Figure 1: Left: The tangent planes at points x and y of the sphere S 2 are different: the vectors v and w of T x M cannot be compared to the vectors t and u of T y M. Thus, it is natural to define the scalar product on each tangent plane. Right: The geodesics starting at x are straight lines in the exponential map and the distance along them is conserved.

Exponential chart

Let x be a point of the manifold that we consider as a local reference and -→ xy a vector of the tangent space T x M at that point. From the theory of second order differential equations, we know that there exists one and only one geodesic starting from that point with this tangent vector. This allows to develop the manifold in the tangent space along the geodesics (think of rolling a sphere along its tangent plane at a given point). The geodesics going through the reference point are transformed into straight lines and the distance along these geodesics is conserved (at least in a neighborhood of x).

The function that maps to each vector -→ xy ∈ T x M the point y of the manifold that is reached after a unit time by the geodesic starting at x with this tangent vector is called the exponential map. This map is defined in the whole tangent space T x M (since the manifold is geodesically complete) but it is generally one-to-one only locally around 0 in the tangent space (i.e. around x in the manifold). In the sequel, we denote by -→ xy = log x (y) the inverse of the exponential map: this is the smallest vector such that y = exp x ( -→ xy). If we look for the maximal definition domain, we find out that it is a star-shaped domain delimited by a continuous curve C x called the tangential cut-locus. The image of C x by the exponential map is the cut locus C x of point x. This is the closure of the set of points where several minimizing geodesics starting from x meet. On the sphere S 2 (1) for instance, the cut locus of a point x is its antipodal point and the tangential cut locus is the circle of radius π.

The exponential map within this domain realizes a chart called the exponential chart. It covers all the manifold except the cut locus of the reference point x, which has a null measure. In this chart, geodesics starting from x are straight lines, and the distance from the reference point are conserved. This chart is somehow the "most linear" chart of the manifold with respect to the reference point x.

Practical implementation

In fact, most of the usual operations using additions and subtractions may be reinterpreted in a Riemannian framework using the notion of bipoint, an antecedent of vector introduced during the 19th Century. Indeed, one defines vectors as equivalent classes of bipoints (oriented couples of points) in a Euclidean space. This is possible because we have a canonical way (the translation) to compare what happens at two different points. In a Riemannian manifold, we can still compare things locally (by parallel transportation), but not any more globally. This means that each "vector" has to remember at which point of the manifold it is attached, which comes back to a bipoint.

However, one can also see a vector -→ xy (attached at point x) as a vector of the tangent space at that point. Such a vector may be identified to a point on the manifold using the geodesic starting at x with tangent vector -→ xy, i.e. using the exponential map: y = exp x ( -→ xy). Conversely, the logarithmic map may be used to map almost any bipoint (x, y) into a vector -→ xy = log x (y) of T x M. This reinterpretation of addition and subtraction using logarithmic and exponential maps is very powerful to generalize algorithms working on vector spaces to algorithms on Riemannian manifolds, as illustrated by Table 1. It is also very powerful in terms of implementation since we can practically express all the geometric operations in these terms: the implementation of log x and exp x is the basis of any programming on Riemannian manifolds, as we will see in the following.

Basic statistical tools

The Riemannian metric induces an infinitesimal volume element on each tangent space, and thus a measure dM on the manifold that can be used to measure random events on the manifold and to define the probability density function (if it exists) of these random elements. It is worth noticing that the induced measure dM represents the notion of uniformity according to the chosen Riemannian metric. This automatic derivation of the uniform measure from the metric gives a rather elegant solution to the Bertrand paradox for geometric probabilities [Poincaré, 1912, Kendall and[START_REF] Kendall | [END_REF]. However, the problem is only shifted: which Riemannian metric do we have to choose? We address this question in Section 3 for real positive definite symmetric matrices (tensors): it turns out that requiring an invariance by the full linear group will lead to a very regular and convenient manifold structure.

Let us come back to the basic statistical tools. With the probability measure of a random element, we can integrate functions from the manifold to any vector space, thus defining the expected value of this function. However, we generally cannot integrate manifold-valued functions. Thus, one cannot define the mean or expected "value" of a random manifold element using a weighted sum or an integral as usual. One solution is to rely on a distancebased variational formulation: the Fréchet or Karcher expected features basically minimize globally (or locally) the variance. As the mean is now defined through a minimization procedure, its existence and uniqueness are not ensured any more (except for distributions with a sufficiently small compact support). In practice, one mean value almost always exists, and it is unique as soon as the distribution is sufficiently peaked. The properties of the mean are very similar to those of the modes (that can be defined as central Karcher values of order 0) in the vectorial case.

To compute the mean value, we designed in [Pennec, 1999, Pennec, 2004] an original Gauss-Newton gradient descent algorithm that essentially alternates the computation of the barycenter in the exponential chart centered at the current estimation of the mean value, and a re-centering step of the chart at the point of the manifold that corresponds to the computed barycenter (geodesic marching step). To define higher moments of the distribution, we used the exponential chart at the mean point: the random feature is thus represented as a random vector with null mean in a star-shaped domain. With this representation, there is no difficulty to define the covariance matrix and potentially higher order moments.

Vector space

Riemannian manifold Subtraction -→

xy = y -x -→ xy = log x (y) Addition y = x + -→ xy y = exp x ( -→ xy) Distance dist(x, y) = y -x dist(x, y) = -→ xy x Mean value (implicit) i -→ xx i = 0 i log x(x i ) = 0 Gradient descent x t+ε = x t -ε∇C(x t ) x t+ε = exp xt (-ε∇C(x t )) Linear (geodesic) interpolation x(t) = x 1 + t --→ x 1 x 2 x(t) = exp x 1 (t --→ x 1 x 2 )
Table 1: Re-interpretation of basic standard operations in a Riemannian manifold.

Based on this covariance matrix, we defined a Mahalanobis distance between a random and a deterministic feature that basically weights the distance between the deterministic feature and the mean feature using the inverse of the covariance matrix. Interestingly, the expected Mahalanobis distance of a random element with itself is independent of the distribution and is equal to the dimension of the manifold, as in the vectorial case.

As for the mean, we chose in [Pennec, 1996, Pennec, 1999, Pennec, 2004] a variational approach to generalize the Normal Law: we define it as the distribution that minimizes the information knowing the mean and the covariance. This amounts to consider a Gaussian distribution on the exponential chart centered at the mean point that is truncated at the cut locus (if there is one). However, the relation between the concentration matrix (the "metric" used in the exponential of the probability density function) and the covariance matrix is slightly more complex than the simple inversion of the vectorial case, as it has to be corrected for the curvature of the manifold. Last but not least, using the Mahalanobis distance of a normally distributed random feature, we can generalize the χ 2 law: we were able to show that is has the same density as in the vectorial case up to an order 3 in σ. This opens the way to the generalization of many other statistical tests, as we may expect similarly simple approximations for sufficiently centered distributions.

Working on the Tensor space

Let us now focus on the space Sym + n of positive definite symmetric matrices (tensors). The goal is to find a Riemannian metric with interesting enough properties. It turns out that it is possible to require an invariance by the full linear group (Section 3.3). This leads to a very regular manifold structure where tensors with null and infinite eigenvalues are both at an infinite distance of any positive definite symmetric matrix: the cone of positive definite symmetric matrices is replaced by a space which has an infinite development in each of its n(n + 1)/2 directions. Moreover, there is one and only one geodesic joining any two tensors, and we can even define globally consistent orthonormal coordinate systems of tangent spaces. Thus, the structure we obtain is very close to a vector space, except that the space is curved.

Exponential, logarithm and square root of tensors

In the following, we will make an extensive use of a few functions on symmetric matrices. The exponential of any matrix can be defined using the series exp(A) = +∞ k=0 A k k! . In the case of symmetric matrices, we have some important simplifications. Let W = U D U T be a diagonalization, where U is an orthonormal matrix, and D = DIAG(d i ) is the diagonal matrix of the eigenvalues. We can write any power of W in the same basis: W k = U D k U T . This means that we may factor out the rotation matrices in the series and map the exponential individually to each eigenvalue:

exp(W ) = +∞ k=0 W k k! = U DIAG(exp(d i )) U T .
The series defining the exponential function converges for any (symmetric) matrix argument, but this is generally not the case for the series defining its inverse function: the logarithm. However, any tensor can be diagonalized into Σ = U DIAG(d i ) U T with strictly positive eigenvalues d i . Thus, the function log(Σ) = U (DIAG(log(d i ))) U T is always well defined on tensors. Moreover, if all the eigenvalues are small enough (|d i -1| < 1), then the series defining the usual log converges and we have:

log(Σ) = U DIAG +∞ k=1 (-1) k+1 k (d i -1) k U T = +∞ k=1 (-1) k+1 k (Σ -Id) k . (1)
The logarithm we defined is obviously the inverse function of exp. Thus, the matrix exponential realizes a one-to-one mapping between the space of symmetric matrices to the space of tensors. Classically, one defines the (left) square root of a matrix B as the set {B

1/2 L } = {A ∈ GL n / A A T = B}. One could also define the right square root: {B 1/2 R } = {A ∈ GL n / A T A = B}.
For tensors, we define the square root as:

Σ 1/2 = {Λ ∈ Sym + n / Λ 2 = Σ}.
The square root is always defined and moreover unique: let Σ = U D 2 U T be a diagonalization (with positives values for the d i 's). Then Λ = U D U T is of course a square root of Σ, which proves the existence. For the uniqueness, let us consider two symmetric and positive square roots Λ 1 and Λ 2 of Σ. Then, Λ 2 1 = Σ and Λ 2 2 = Σ obviously commute and thus they can be diagonalized in the same basis: this means that the diagonal matrices D 2 1 and D 2 2 are equal. As the elements of D 1 and D 2 are positive, they are also equal and Λ 1 = Λ 2 . Last but not least, we have the property that Σ 1/2 = exp 1 2 (log Σ) .

An affine invariant distance

Let us consider the following action of the linear group GL n on the tensor space Sym + n :

A Σ = AΣA T ∀A ∈ GL n and Σ ∈ Sym + n .
This group action corresponds for instance to the standard action of the affine group on the covariance matrix Σ xx of a random variables

x in R n : if y = Ax + t, then ȳ = Ax + t and Σ yy = E[(y -ȳ) (y -ȳ) T ] = AΣ xx A T .
This action is naturally extended to tangent vectors is the same way

: if Γ(t) = Σ + t W + O(t 2 ) is a curve passing at Σ with tangent vector W , then the curve A Γ(t) = A Σ A T + t A W A T + O(t 2 ) passes through A Σ with tangent vector A W .
Following [START_REF] Pennec | Uniform distribution, distance and expectation problems for geometric features processing[END_REF]], any invariant distance on

Sym + n verifies dist(A Σ 1 , A Σ 2 ) = dist(Σ 1 , Σ 2 ). Choosing A = Σ -1/2 1
, we can reduce this to a pseudo-norm, or distance to the identity:

dist(Σ 1 , Σ 2 ) = dist Id, Σ -1 2 1 Σ 2 Σ -1 2 1 = N Σ -1 2 1 Σ 2 Σ -1 2 1 .
Moreover, as the invariance has to hold for any transformation, N should be invariant under the action of the isotropy group

H( Id) = O n = {U ∈ GL n / U U T = Id}: ∀U ∈ O n , N (U Σ U T ) = N (Σ).
Using the spectral decomposition Σ = U D 2 U T , it is easy to see that N (Σ) has to be a symmetric function of the eigenvalues. Moreover, the symmetry of the distance dist(Σ, Id) = dist( Id, Σ) imposes that N (Σ) = N (Σ (-1) ). Thus, a good candidate is the sum of the squared logarithms of the eigenvalues:

N (Σ) 2 = log(Σ) 2 = n i=1 (log(σ i )) 2 .
(2)

This "norm" verifies by construction the symmetry and positiveness. N (Σ) = 0 implies that σ i = 1 (and conversely), so that the separation axiom is verified. However, we do no know any simple proof of the triangle inequality, which should read

N (Σ 1 ) + N (Σ 2 ) ≥ N (Σ -1/2 1 Σ 2 Σ -1/2 1
), even if we can verify it experimentally (see e.g. [START_REF] Förstner | A metric for covariance matrices[END_REF]).

An invariant Riemannian metric

Another way to determine the invariant distance is through the Riemannian metric. Let us take the most simple scalar product on the tangent space at the identity matrix: if W 1 and W 2 are tangent vectors (i.e. symmetric matrices, not necessarily definite nor positive), we define the scalar product to be the standard matrix scalar product

W 1 | W 2 = Tr(W T 1 W 2
). This scalar product if obviously invariant by the isotropy group O n . Now, if W 1 and W 2 are two tangent vectors at Σ, we require their scalar product to be invariant by the action of any transformation:

W 1 | W 2 Σ = A W 1 | A W 2 A Σ .
This should be true in particular for A = Σ -1/2 , which allows us to define the scalar product at any Σ from the scalar product at the identity:

W 1 | W 2 Σ = Σ -1 2 W 1 Σ -1 2 Σ -1 2 W 2 Σ -1 2 Id = Tr Σ -1 2 W 1 Σ -1 W 2 Σ -1 2 .
One can easily verify that this definition is left unchanged if we use any other transformation A = U Σ -1/2 (where U is a free orthonormal matrix) that transports Σ to the identity:

A Σ = A Σ A T = U U T = Id.
To find the geodesic without going though the computation of Christoffel symbols, we may rely on a result from differential geometry [Gamkrelidze, 1991, Helgason, 1978, Kobayashi and Nomizu, 1969 which says that the geodesics for the invariant metrics on affine symmetric spaces are generated by the action of the one-parameter subgroups of the acting Lie group 1 . Since the one-parameter subgroups of the linear group are given by the matrix exponential exp(t A), geodesics on our tensor manifold going through Σ with tangent vector W should have the following form:

Γ (Σ,W ) (t) = exp(t A) Σ exp(t A) T with W = A Σ + Σ A T . (3)
For our purpose, we need to relate explicitly the geodesic to the tangent vector in order to define the exponential chart. Since Σ is a symmetric matrix, there is hopefully an explicit solution to the Sylvester equation

W = A Σ + Σ A T . We get A = 1 2 W Σ (-1) + Σ 1/2 Z Σ -1/2
, where Z is a free skew-symmetric matrix. However, introducing this solution into the equation of geodesics (Eq. 3) does not lead to a very tractable expression. Let us look at an alternative solution.

Since our metric (and thus the geodesics) is invariant under the action of the group, we can focus on the geodesics going through the origin (the identity). In that case, a symmetric solution of the Sylvester equation is A = 1 2 W , which gives the following equation for the geodesic going through the identity with tangent vector W :

Γ ( Id,W ) (t) = exp t 2 W exp t 2 W T = exp(t W ).
We may observe that the tangent vector along this curve is the parallel transportation of the initial tangent vector.

If W = U DIAG(w i ) U T , dΓ(t) dt = U DIAG (w i exp(t w i )) U T = Γ(t) 1 2 W Γ(t) 1 2 = Γ(t) 1 2 W.
By definition of our invariant metric, the norm of this vector is constant:

Γ(t) 1/2 W 2 Γ(t) 1/2 Id = W 2 Id = W 2 2 .
This was expected since geodesics are parameterized by arc-length. Thus, the length of the curve between time 0 and 1 is

L = 1 0 dΓ(t) dt 2 Γ(t) dt = W 2 Id .
Solving for Γ ( Id,W ) (1) = Σ, we obtain the "norm" N (Σ) of Eq.( 2). Using the invariance of our metric, we easily obtain the geodesic starting from any other point of the manifold using our group action:

Γ (Σ,W ) (t) = Σ 1 2 Γ ( Id,Σ -1/2 W ) (t) = Σ 1 2 exp t Σ -1 2 W Σ -1 2 Σ 1 2 .
Coming back to the distance dist 2 (Σ, Id) = i (log σ i ) 2 , it is worth noticing that tensors with null eigenvalues are located as far from the identity as tensors with infinite eigenvalues: at the infinity. Thanks to the invariance by the linear group, this property holds for the distance to any (positive definite) tensor of the manifold. Thus, the original cone of positive definite symmetric matrices (a linear manifold with a flat metric but which is incomplete: there is a boundary at a finite distance) has been changed into a regular and complete (but curved) manifold with an infinite development in each of its n(n + 1)/2 directions.

Exponential and logarithm maps

As a general property of Riemannian manifolds, geodesics realize a local diffeomorphism from the tangent space at a given point of the manifold to the manifold: Γ (Σ,W ) (1) = exp Σ (W ) associates to each tangent vector W ∈ T Σ Sym + n a point of the manifold. This mapping is called the exponential map, because it corresponds to the usual exponential in some matrix groups. This is exactly our case for the exponential map around the identity:

exp Id (U DU T ) = exp(U DU T ) = U DIAG (exp(d i )) U T .
However, the Riemannian exponential map associated to our invariant metric has a more complex expression at other tensors:

exp Σ (W ) = Σ 1 2 exp Σ -1 2 W Σ -1 2 Σ 1 2 .
In our case, this diffeomorphism is global, and we can uniquely define the inverse mapping everywhere:

log Σ (Λ) = Σ 1 2 log Σ -1 2 ΛΣ -1 2 Σ 1 2 .
Thus, exp Σ gives us a collection of one-to-one and complete maps of the manifold, centered at any point Σ. As explained in Section 2.1, these charts can be viewed as the development of the manifold onto the tangent space along the geodesics. Moreover, as the manifold has a non-positive curvature [Skovgaard, 1984], there is no cut-locus and the statistical properties detailed in [Pennec, 2004] hold in their most general form. For instance, we have the existence and uniqueness of the mean of any distribution with a compact support [Kendall, 1990].

Induced and orthonormal coordinate systems

One has to be careful because the coordinate system of all these charts is not orthonormal. Indeed, the coordinate system of each chart is induced by the standard coordinate system (here the matrix coefficients), so that the vector -→ ΣΛ corresponds to the standard derivative in the vector space of matrices: we have Λ = Σ + -→ ΣΛ + O( -→ ΣΛ2 ). Even if this basis is orthonormal at some points of the manifold (such as at the identity for our tensors), it has to be corrected for the Riemannian metric at other places due to the manifold curvature.

From the expression of the metric, one can observe that

-→ ΣΛ 2 Σ = log Σ (Λ) 2 Σ = Σ -1 2 log Σ (Λ)Σ -1 2 2 Id = log(Σ -1 2 Λ) 2 2 .
This shows that

-→ ΣΛ ⊥ = log(Σ -1 2 Λ) ∈ T Σ Sym + n is the expression of the vector -→
ΣΛ in an orthonormal basis. In our case, the transformation Σ 1/2 ∈ GL n is moreover uniquely defined (as a positive square root) and is a smooth function of Σ over the complete tensor manifold. Thus, -→ ΣΛ ⊥ realizes an atlas of orthonormal exponential charts which is globally smooth with respect to the reference point 2 Σ. This group action approach was chosen in earlier works [Pennec, 1996, Pennec and Thirion, 1997, Pennec and Ayache, 1998] with what we called the placement function.

For some statistical operations, we need to use a minimal representation (e.g. 6 parameters for 3 × 3 tensors) in a (locally) orthonormal basis. This can be realized through the classical "Vec" operator that maps the element a i,j of a n × n matrix A to the (i n + j)th element Vec(A) i n+j of a n × n dimensional vector Vec(A). Since we are working with symmetric matrices, we have only n(n + 1)/2 independent coefficients, say the upper triangular part. However, the off-diagonal coefficients are counted twice in the L 2 norm at the identity: W 2 2 = n i=1 w 2 i,i + 2 i<j≤n w 2 i,j . Thus, to express our minimal representation in an orthonormal basis, we need to multiply the off diagonal terms by √ 2:

Vec Id (W ) = w 1,1 , √ 2 w 1,2 , w 2,2 , √ 2 w 1,3 , √ 2 w 2,3 , w 3,3 , . . . √ 2 w 1,n , . . . √ 2 w (n-1),n , w n,n T .
Now, for a vector -→ ΣΛ ∈ T Σ Sym + n , we define its minimal representation in the orthonormal coordinate system as:

Vec Σ ( -→ ΣΛ) = Vec Id ( -→ ΣΛ ⊥ ) = Vec Id Σ -1 2 -→ ΣΛ Σ -1 2 = Vec Id log(Σ -1 2 Λ) .
The mapping Vec Σ realizes an explicit isomorphism between T Σ Sym + n and R n(n+1)/2 with the canonical metric.

Gradient descent and PDEs: an intrinsic geodesic marching scheme

Let f (Σ) be an objective function to minimize, Σ t the current estimation of Σ, and

W t = ∂ Σ f = [∂f /∂σ ij ]
its matrix derivative at that point, which is of course symmetric. The principle of a first order gradient descent is to go toward the steepest descent, in the direction opposite to the gradient for a short time-step ε, and iterate the process. However, the standard operator Σ t+1 = Σ t -εW t is only valid for very short time-steps in the flat Euclidean matrix space, and we could easily go out of the cone of positive definite tensors. A much more interesting numerical operator is given by following the geodesic backward starting at Σ with tangent vector W t during a time ε. This intrinsic gradient descent ensures that we cannot leave the manifold. It can easily be expressed using the exponential map:

Σ t+1 = Γ (Σt,Wt) (-ε) = exp Σt (-εW t ) = Σ 1 2 exp(-εΣ -1 2 W t Σ -1 2 )Σ 1 2 .
This intrinsic scheme is trivially generalized to partial differential evolution equations (PDEs) on tensor fields such as

∂ t Σ(x, t) = -W (x, t): we obtain Σ(x, t+dt) = exp Σ(x,t) (-dtW (x, t)).

Example with the mean value

Let Σ 1 . . . Σ N be a set of measures of the same Tensor. The Karcher or Fréchet mean is the set of tensors minimizing the sum of squared distances:

C(Σ) = N i=1 dist 2 (Σ, Σ i ).
In the case of tensors, the manifold has a non-positive curvature [Skovgaard, 1984], so that there is one and only one mean value Σ [Kendall, 1990]. Moreover, a necessary and sufficient condition for an optimum is a null gradient of the criterion. Thus, the intrinsic Newton gradient descent algorithm gives the following mean value at estimation step t + 1:

Σt+1 = expΣ t 1 N N i=1 logΣ t (Σ i ) = Σ 1 2 t exp 1 N N i=1 log Σ-1 2 t Σ i Σ-1 2 t Σ 1 2 t . (4)
Note that we cannot easily simplify this expression further as in general the data Σ i and the mean value Σt cannot be diagonalized in a common basis. However, this gradient descent algorithm usually converges very fast (about 10 iterations, see Fig. 2).

Simple statistical operations on tensors

As described in [Pennec, 2004], we may generalize most of the usual statistical methods by using the exponential chart at the mean point. For instance, the empirical covariance matrix of a set of N tensors Σ i of mean Σ is defined using the tensor product:

1 N -1 n i=1 --→ ΣΣ i ⊗ --→ ΣΣ i .
Using our Vec mapping, we may come back to more usual matrix notations and write its expression in a minimal representation with an orthonormal coordinate system:

Cov = 1 N -1 N i=1 VecΣ --→ ΣΣ i VecΣ --→ ΣΣ i T .
One may also define the Mahalanobis distance

µ 2 ( Σ,Cov) (Σ) = VecΣ -→ ΣΣ T Cov (-1) VecΣ -→ ΣΣ .
Looking for the probability density function that minimizes the information with a constrained mean and covariance, we obtain a generalization of the Gaussian distribution of the form:

NΣ ,Γ (Σ) = k exp - 1 2 µ 2 Σ,Γ (Σ) .
The main difference with a Euclidean space is that we have a curvature to take into account: the invariant measure induced on the manifold by our metric is linked to the usual matrix measure by dM(Σ) = dΣ/ det(Σ). Likewise, the curvature slightly modifies the usual relation between the covariance matrix, the concentration matrix Γ and the normalization parameter k of the Gaussian distribution [Pennec, 2004]. These differences have an impact on the calculations using continuous probability density functions. However, from a practical point of view, we only deal with a discrete sample set of measurements, so that the measure-induced corrections are hidden. For instance, we can generate a random (generalized) Gaussian tensor using the following procedure: we sample n(n + 1)/2 independent and normalized real Gaussian samples, multiply the corresponding vector by the square root of the desired covariance matrix (expressed in our Vec coordinate system), and come back to the tensor manifold using the inverse Vec mapping. Using this procedure, we can easily generate noisy measurements of known tensors (see e.g. Fig. 7).

To check the implementation of our charts and geodesic marching algorithms, we verified experimentally the central limit theorem. This theorem states that the empirical mean of N independently and identically distributed (IID) random variables with a variance γ 2 asymptotically follows a Gaussian law of variance γ 2 /N , centered at the exact mean value. The principle of our experiments is now as follows. We randomly generated N random Gaussian tensors around a random tensor Σ with a variance of γ 2 = 1. We computed the mean Σ using the algorithm of Eq. 4. The convergence is clearly very fast (Fig. 2, left). Now, if the error between the computed and the exact mean really follows a Gaussian law of variance γ 2 /N , then the normalized Mahalanobis distance µ 2 = N dist( Σ, Σ) 2 /γ 2 should follow a χ 2 6 distribution. However, this simple experiment only gives us one measurement. Thus, to verify the distribution, we repeated this experiment with N varying from 10 to 1000. Figure 2 presents the histogram of the normalized Mahalanobis distances we obtain. The empirical distribution follows quite well the theoretical χ 2 6 distribution, as expected, with an empirical mean of 5.85 and a variance of 12.17 (expected values are 6 and 12). Moreover, a Kolmogorov-Smirnov test confirms that the distance between the empirical and theoretical cumulative pdf is not significant (p-value of 0.19).

Tensor Interpolation

One of the important operations in geometric data processing is to interpolate values between known measurements. In 3D image processing, (tri-) linear interpolation is often used thanks to its very low computational load and comparatively much better results than nearest neighbor interpolation. Other popular methods include the cubic and, more generally, spline interpolations [START_REF] Thévenaz | Interpolation revisited[END_REF], Meijering, 2002].

The standard way to define an interpolation on a regular lattice of dimension d is to consider that the interpolated function f (x) is a linear combination of samples

f k at integer (lattice) coordinates k ∈ Z d : f (x) = k w(x -k) f k .
To realize an interpolation, the "sample weight" function w has to vanish at all integer coordinates except 0 where it has to be one. A typical example where the convolution kernel has an infinite support is the sinus cardinal interpolation. With the nearest-neighbor, linear (or tri-linear in 3D), and higher order spline interpolations, the kernel is piecewise polynomial, and limited to a few neighboring points in the lattice.

When it comes to an irregular sampling (i.e. a set of measurements f k at positions x k ), interpolation may still be defined using a weighted mean:

f (x) = N k=1 w k (x) f k .
To ensure that this is an interpolating function, one has to require that w i (x j ) = δ ij (where δ ij is the Kronecker symbol). Moreover, the coordinates are usually normalized so that N k=1 w k (x) = 1 for all position x within the domain of interest. Typical examples in triangulations or tetrahedrizations are barycentric and natural neighbor coordinates [Sibson, 1981].

Interpolation through weighted mean

To generalize interpolation methods defined using weighted means to our tensor manifold, let us assume that the sample weights w k (x) are defined as above in R d . Thanks to their normalization, the value f (x) interpolated from vectors f k verifies N i=1 w i (x)(f i -f (x)) = 0. Thus, similarly to the Fréchet mean, we can define the interpolated value Σ(x) on our tensor manifold as the tensor that minimizes the weighted sum of squared distances to the measurements Σ i :

C(Σ(x)) = N i=1 w i (x) dist 2 (Σ i , Σ(x)).
Of course, we loose in general the existence and uniqueness properties. However, for positive weights, the existence and uniqueness theorems for the Karcher mean can be adapted. In practice, this means that we have a unique tensor that verifies N i=1 w i (x)

----→ Σ(x)Σ i = 0. To reach this solution, it is easy to adapt the Gauss-Newton scheme proposed for the Karcher mean. The algorithm becomes:

Σ t+1 (x) = exp Σt(x) N i=1 w i (x) log Σt(x) (Σ i ) = Σ 1 2 t (x) exp N i=1 w i (x) log Σ -1 2 t (x)Σ i Σ -1 2 t (x) Σ 1 2 t (x). (5)
Once again, this expression cannot be easily simplified, but the convergence is very fast (usually less than 10 iterations as for the mean).

Example of the linear interpolation

The linear interpolation is simple as this is a walk along the geodesic joining the two tensors.

For instance, the interpolation in the standard Euclidean matrix space would give Σ (t) = (1 -t) Σ 1 + t Σ 2 . In our Riemannian space, we have the closed-form expression:

Σ(t) = exp Σ 1 (t log Σ 1 (Σ 2 )) = exp Σ 2 ((1 -t) log Σ 2 (Σ 1 )) for t ∈ [0; 1].
We displayed in Fig. 3 the flat and the Riemannian interpolations between 2D tensors of eigenvalues (5,1) horizontally and (1,50) at 45 degrees, along with the evolution of the eigenvalues, their mean (i.e. trace of the matrix) and product (i.e. determinant of the matrix or volume of the ellipsoid).

With the standard matrix coefficient interpolation, the evolution of the trace is perfectly linear (which was expected since this is a linear function of the coefficients), and the principal eigenvalue regularly grows almost linearly, while the smallest eigenvalue slightly grows toward a local maxima before lowering. What is much more annoying is that the determinant (i.e. the volume) does not grow regularly in between the two tensors, but goes through a maximum. If we interpret our tensors as covariance matrices of Gaussian distributions, this means that the probability of a random point to be accepted as a realization of our distribution is larger in between than at the measurement points themselves! On the contrary, one can clearly see a regular evolution of the eigenvalues and of their product with the interpolation in our Riemannian space. Moreover, there is a much smoother rotation of the eigenvectors than with the standard interpolation.

Tri-linear interpolation

The bi-and tri-linear interpolation of tensors on a regular grid in 2D or 3D are almost as simple, except that we do not have any longer an explicit solution using geodesics since there are more than two reference points. After computing the (bi-) tri-linear weights with respect to the neighboring sites of the point we want to evaluate, we now have to go through the iterative optimization of the weighted mean (Eq. 5) to compute the interpolated tensor. We display an example in Figure 4. One can see that the volume of the tensors is much more important with the classical than with the Riemannian interpolation. We also get a much smoother interpolation of the principal directions with our method.

Interpolation of non regular measurements

When tensors are not measured on a regular grid but "randomly" localized in space, defining neighbors becomes an issue. One solution, proposed by [Sibson, 1981] and later used for surfaces by [START_REF] Cazals | Natural coordinates of points on a surface[END_REF], is the natural neighbor interpolation. For any point x, its natural neighbors are the points of {x i } whose Voronoi cells are chopped off upon insertion of x into the Voronoi diagram. The weight w i of each natural neighbor x i is the proportion of the new cell that is taken away by x to x i in the new Voronoï diagram. One important restriction of these interesting coordinates is that they are limited to the convex hull of the point set (otherwise the volume or surface of the cell is infinite).

Another idea is to rely on radial-basis functions to define the relative influence of each measurement point. For instance, a Gaussian influence would give a weight w i (x) = G σ (xx i ) to the measurement Σ i located at x i . Since weights need to be renormalized in our setup, this would lead to the following evolution equation:

Σ t+1 (x) = exp Σt(x) N i=1 G σ (x -x i ) -----→ Σ t (x)Σ i N i=1 G σ (x -x i ) . ( 6
)
The initialization could be the (normalized) Gaussian mean in the matrix space. An example of the result of this evolution scheme is provided on top of Figure 10. However, this algorithm does not lead to an interpolation, but rather to an approximation, since the weights are not zero at other measurement points. Moreover, we have little control on the quality of this approximation. It is only at the limit where σ goes to zero that we end-up with a (noncontinuous) closest point interpolation. We will describe in Section 6.3 a last alternative that performs the interpolation and extrapolation of sparsely distributed tensor measurements using diffusion.

Filtering tensor fields

Let us now consider that we have a tensor field, for instance like in Diffusion Tensor Imaging (DTI) [Le [START_REF] Bihan | [END_REF], where the tensor is a first order approximation of the anisotropic diffusion of the water molecules at each point of the imaged tissues. In the brain, the diffusion is much favored in the direction of oriented structures (fibers of axons). One of the goal of DTI is to retrieve the main tracts along these fibers. However, the tensor field obtained from the images is noisy and needs to be regularized before being further analyzed. A naive but simple and often efficient regularization on signal or images is the convolution by a Gaussian. The generalization to tensor fields is quite straightforward using once again weighted means (Section 5.1 below). An alternative is to consider a regularization using diffusion. This will be the subject of Sections 5.3 and 5.4.

Gaussian Filtering

In the continuous setting, the convolution of a vector field F 0 (x) by a Gaussian is:

F (x) = y G σ (y -x) F 0 (y) dy.
In the discrete setting, coefficients are renormalized since the neighborhood V is usually limited to points within one to three times the standard deviation:

F (x) = u∈V(x) G σ (u) F 0 (x + u) u∈V(x) G σ (u) = arg min F u∈V(x) G σ (u) F 0 (x + u) -F 2 .
Like previously, this weighted mean can be solved on our manifold using our intrinsic gradient descent scheme. Starting from the measured tensor field Σ 0 (x), the evolution equation is

Σ t+1 (x) = exp Σt(x) u∈V G σ (u) ----------→ Σ t (x)Σ t (x + u) u∈V G σ (u)
.

We illustrate in Fig. 5 the comparative Gaussian filtering of a slice of a DT MR image using the flat metric on the coefficient (since weights are positive, a weighted sum of positive definite matrices is still positive definite) and our invariant Riemannian metric. One can see a more important blurring of the corpus callosum fiber tracts using the flat metric.

However, the integration of this filtering scheme into a complete fiber tracking system would be necessary to fully evaluate the pros and cons of each metric. Left: raw estimation of the tensors. The color codes for the direction of the principal eigenvector (red: left/right, green: anterior/posterior, blue: top/bottom). Middle: Gaussian filtering of the coefficients (5x5 window, σ = 2.0). Right: equivalent filtering (same parameters) using the Riemannian metric.

Spatial gradient of Tensor fields

On a n-dimensional vector field F (x) = (f 1 (x 1 , . . . x d ), . . . f n (x 1 , . . . x d )) T over R d , one may express the spatial gradient in an orthonormal basis as:

∇F T = ∂F ∂x = [∂ 1 F, . . . ∂ d F ] =    ∂f 1 ∂x 1 , . . . ∂f 1 ∂x d . . . . . . . . . ∂fn ∂x 1 , . . . ∂fn ∂x d    .
The linearity of the derivatives implies that we could use directional derivatives in more than the d orthogonal directions. This is especially well adapted to stabilize the discrete computations: the finite difference estimation of the directional derivative is

∂ u F (x) = F (x+ u) -F (x)
. By definition, the spatial gradient is related to the directional derivatives through ∇F T u = ∂ u F (x). Thus, we may compute ∇F as the matrix that best approximates (in the least-square sense) the directional derivatives in the neighborhood V (e.g. 6, 18 or 26 connectivity in 3D):

∇F (x) = arg min G u∈V G T u -∂ u F (x) 2 = u∈V u u T (-1) u∈V u ∂ u F (x) T u∈V u u T (-1) u∈V u (F (x + u) -F (x)) T .
We experimentally found in other applications (e.g. to compute the Jacobian of a deformation field in non-rigid registration [Rey et al., 2002, p. 169]) that this gradient approximation scheme was more stable and much faster than computing all derivatives using convolutions, for instance by the derivative of the Gaussian.

To quantify the local amount of variability independently of the space direction, one usually takes the norm of the gradient:

∇F (x) 2 = d i=1 ∂ i F (x) 2 .
Once again, this can be approximated using all directional derivatives in the neighborhood

∇F (x) 2 d Card(V) u∈V F (x + u) -F (x) 2 u 2 . ( 7
)
Notice that this approximation is consistent with the previous one only if the directions u are normalized to unity. For a manifold valued field Σ(x) define on R d , we can proceed similarly, except that the directional derivatives ∂ i Σ(x) are now tangent vectors of T Σ(x) M. They can be approximated just like above using finite "differences" in our exponential chart:

∂ u Σ(x) ---------→ Σ(x) Σ(x + u) = Σ(x) 1 2 log Σ(x) -1 2 Σ(x + u) Σ(x) -1 2 Σ(x) 1 2 . ( 8
)
As observed in Section 3.5, we must be careful that this directional derivative is expressed in the standard matrix coordinate system (coefficients). Thus, the basis is not orthonormal: to quantify the local amount of variation, we have to take the metric at the point Σ(x) into account, so that:

∇Σ(x) 2 Σ(x) = d i=1 ∂ i Σ(x) 2 Σ(x) d Card(V) u∈V log Σ(x) -1 2 Σ(x + u) Σ(x) -1 2 2 2 u 2 . (9)

Filtering using PDEs

Regularizing a scalar, vector or tensor field F aims at reducing the amount of its spatial variations. The first order measure of such variations is the spatial gradient ∇F that we dealt with in the previous section. To obtain a regularity criterion over the domain Ω, we just have to integrate:

Reg(F ) = Ω ∇F (x) 2 dx.
Starting from an initial field F 0 (x), the goal is to find at each step a field F t (x) that minimizes the regularity criterion by gradient descent in the space of (sufficiently smooth and square integrable) functions.

To compute the first order variation, we write a Taylor expansion for an incremental step in the direction of the field H. Notice that H(x) is a tangent vector at F (x):

Reg(F + ε H) = Reg(F ) + 2 ε Ω ∇F (x) | ∇H(x) dx + O(ε 2 ).
We get the directional derivative: ∂ H Reg(F ) = 2 Ω ∇F (x) | ∇H(x) dx. To compute the steepest descent, we now have to find the gradient ∇Reg(F ) such that for all variation H, we have ∂ H Reg(F ) = Ω ∇Reg(F )(x) | H(x) F (x) dx. Notice that ∇Reg(F )(x) and H(x) are elements of the tangent space at F (x), so that the scalar product should be taken at F (x) for a tensor field. (with the flat metric). Middle: we computed the directional derivatives with the exponential map (Eq. 8), but the norm is taken without correcting for the metric. As this should be very close to the flat gradient norm, we only display the difference image. The main differences are located on very sharp boundaries, where the curvature of our metric has the most important impact. However, the relative differences remains small (less than 10%), which shows the stability of both the gradient and the log / exp computation schemes. Right: Riemannian norm of the Riemannian gradient (Eq. 9). One can see much more detailed structures within the brain, which will now be preserved during an anisotropic regularization step.

The case of a scalar field Let f : R d → R be a scalar field. Our regularization criterion is Reg(f ) = Ω ∇f (x) 2 dx. Let us introduce the divergence div(.) = ∇ | . and the Laplacian operator ∆f = div(∇f ). The divergence is usually written ∇ T = (∂ 1 , . . . , ∂ d ), so that in an orthonormal coordinate system we have ∆f = ∇ | ∇f = d i=1 ∂ 2 i f . Let now G(x) be a vector field. Typically, we will use G(x) = ∇f (x). Using the standard differentiation rules, we have:

div(h G) = ∇ | h G = h div(G) + ∇h | G .
Now, thanks to the Green's formula (see e.g. [START_REF] Gallot | Riemannian Geometry[END_REF]), we know that the flux going out of the boundaries of a (sufficiently smooth) region Ω is equal to the integral of the divergence inside this region. If we denote by n the normal pointing outward at a boundary point, we have:

∂Ω h G | n dn = Ω div(h G) = Ω h div(G) + Ω ∇h | G .
This result can also be interpreted as an integration by part in R d . Assuming homogeneous Neumann boundary conditions (gradient orthogonal to the normal on ∂Ω: G | n = 0), the flow across the boundary vanishes, and we are left with: Ω G | ∇h = -Ω h div(G). Thus, coming back to our original problem, we have:

∂ h Reg(f )(x) = 2 Ω ∇f (x) | ∇h(x) dx = -2 Ω h(x) ∆f (x) dx.
Since this last formula is no more than the scalar product on the space L 2 (Ω, R) of square integrable functions, we end-up with the classical Euler-Lagrange equation: ∇Reg(f ) = -2∆f (x). The evolution equation used to filter the data is thus

f t+1 (x) = f t (x) -ε∇Reg(f )(x) = f t (x) + 2 ε∆f t (x).
The vector case Let us decompose our vector field F (x) into its n scalar components f i (x). Likewise, we can decompose the d × n gradient ∇F into the gradient of the n scalar components ∇f i (x) (columns). Thus, choosing an orthonormal coordinate system on the space R n , our regularization criterion is decoupled into n independent scalar regularization problems:

Reg(F )(x) = n i=1 Ω ∇f i (x) 2 dx = n i=1 Reg(f i ).
Thus, each component f i has to be independently regularized with the Euler-Lagrange equation: ∇Reg(f i ) = -2∆f i . With the convention that the Laplacian is applied component-wise (so that we still have ∆F = div(∇F ) = ∇ T ∇F = (∆f 1 , . . . ∆f n ) T ), we end-up with the vectorial equation:

∇Reg(F ) = -2∆F for Reg(F ) = Ω ∇F (x) dx.
The associated evolution equation is

F t+1 (x) = F t (x) + 2 ε∆F t (x).
Tensor fields For a tensor field Σ(x) ∈ Sym + n over R d , the procedure is more complex as we should use the covariant derivative (the connection) to differentiate vectors fields on our manifold. However, we may avoid the introduction of additional complex mathematical tools by coming back to the basic definitions. We summarize below the main ideas, while the full calculations are worked out in Appendix A.1. Let (x 1 , . . . x d ) be an orthonormal coordinate system of R d . Our regularization criterion is:

Reg(Σ) = Ω ∇Σ(x) 2 Σ(x) dx = d i=1 Ω ∂ i Σ 2 Σ . ( 10
)
The idea is to write this criterion as the trace of sums and products of standard Euclidean matrices and to compute its directional derivative ∂ W Reg for a perturbation field W . This expression contains of course derivatives ∂ i W that we need to integrate. However, as everything is expressed in the standard Euclidean chart (matrix coefficients), and assuming the proper Neumann boundary conditions, we shall safely use the previous integration by part formula Ω Tr((

∂ i W ) Λ i ) = -Ω Tr(W (∂ i Λ i ))
. Notice that we are using the matrix coefficients only as a chart and not as a metric. Eventually, we rewrite the obtained expression in terms of our Riemannian metric to obtain the formula defining the gradient of the criterion:

∂ W Reg = Ω W | ∇Reg Σ .
By identification, we get: ∇Reg(Σ) = -2 ∆Σ, where ∆ is the Laplace-Beltrami operator on our manifold:

∆Σ = d i=1 ∆ i Σ with ∆ i Σ = ∂ 2 i Σ -(∂ i Σ) Σ (-1) (∂ i Σ). (11)
As we can see, the flat Euclidean second order directional derivatives ∂ 2 i Σ are corrected by an additional term due to the curvature of our manifold. To conclude, the gradient descent on the regularization criterion with the intrinsic geodesic marching scheme of Section 3.6 leads to:

Σ t+1 (x) = exp Σt(x) (-ε ∇Reg(Σ)(x)) = exp Σt(x) (2 ε ∆Σ(x)) .
For the numerical computation of the Laplacian, we may approximate the first and second order tensor derivative by their Euclidean derivatives. This gives a fourth order approximation of the Laplace-Beltrami operator (see Appendix A.2). However, this numerical scheme is extrinsic since it is based on (Euclidean) differences of tensors. We propose here an intrinsic scheme based on the exponential chart at the current point. We already know from Eq. 8 that ---------→ Σ(x)Σ(x + u) is an approximation of the first order directional derivative

∂ u Σ(x). We show in Appendix A.2 that ---------→ Σ(x)Σ(x + u) + ---------→ Σ(x)Σ(x -u
) is a forth order approximation of the Laplace Beltrami operator in the direction u:

∆ u Σ = ∂ 2 u Σ -2 (∂ u Σ) Σ (-1) (∂ u Σ) = ---------→ Σ(x)Σ(x + u) + ---------→ Σ(x)Σ(x -u) + O( u 4 ). ( 12
)
To compute the complete manifold Laplacian of Eq. 11, we just have to compute the above numerical approximations of the tensor field derivatives along d orthonormal basis vectors x i . However, like for the computation of the gradient, we may improve the stability of the numerical scheme by averaging the derivatives in all possible directions in the neighborhood V. Assuming a symmetric and isotropic neighborhood, we finally obtain:

∆Σ(x) = d Card(V) u∈V ∆ 2 u Σ(x) u 2 2 d Card(V) u∈V ---------→ Σ(x)Σ(x + u) u 2 . ( 13
)

Anisotropic filtering

In practice, we would like to filter within the homogeneous regions, but not across their boundaries. The basic idea is to penalize the smoothing in the directions where the derivative is important [Perona andMalik, 1990, Gerig et al., 1992]. If c(.) is a weighting function decreasing from c(0) = 1 to c(+∞) = 0, this can be realized directly in the discrete implementation of the Laplacian (Eq. 13): the contribution ∆ u Σ of the spatial direction u to the Laplace-Beltrami operator is weighted by our decreasing function according to the norm ∂ u Σ Σ of the gradient in that direction. The important point here is that we should evaluate the norm of directional derivatives of the tensor field with our invariant metric.

With our finite difference approximations, this leads to the following modified Laplacian:

∆ aniso Σ(x) = d Card(V) u∈V c ∂ u Σ(x) Σ(x) u ∆ 2 u Σ(x) u 2 2 d Card(V) u∈V c    ---------→ Σ(x)Σ(x + u) Σ(x) u    ---------→ Σ(x)Σ(x + u) u 2 . ( 14
)
Figures 7 and 8 present example results of this very simple anisotropic filtering scheme on synthetic and real DTI images. We used the function c(x) = exp (-x 2 /κ 2 ), where the threshold κ controls the amount of local regularization: for a gradient magnitude greater than 2 to 3 times κ, there is virtually no regularization, while the field is almost linearly smoothed for gradient magnitudes below a fraction (say 0.1) of κ. For both synthetic and real data, the histogram of the gradient norm is very clearly bimodal so that the threshold κ is easily determined.

In Fig. 7, we generated a tensor field with a discontinuity, and add independent Gaussian noises according to Section 3.8. The anisotropic smoothing perfectly preserves the discontinuity while completely smoothing each region. In this synthetic experiment, we retrieve tensor values that are very close to the initial tensor field. This could be expected since the two regions are perfectly homogeneous. After enough regularization steps, each region is a constant field equal to the mean of the 48 initially noisy tensors of the region. Thus, similarly to the Euclidean mean of identically and independently distributed measurements, we expect the standard deviation of the regularized tensors to be roughly 7 √ 48 times smaller than the one of the noisy input tensors. In Figure 8, we display the evolution of (a slice of) the tensor field, the norm of the gradient and the fractional anisotropy (FA) at different steps of the anisotropic filtering of a 3D DTI. The FA is based on the normalized variance of the eigenvalues. It shows the differences between an isotropic diffusion in the brain (where the diffusion tensor is represented by a sphere, FA=0) and a highly directional diffusion (cigar-shaped ellipsoid, FA=1). Consequently, the bright regions in the image are the potential areas where nervous fibers are located. One can see that the tensors are regularized in "homogeneous" regions (ventricles, temporal areas), while the main tracts are left unchanged. It is worth noticing that the fractional anisotropy is very well regularized even though this measure has almost nothing in common with our invariant tensor metric.

Figure 9 displays closeups around the ventricles to compare the different regularization methods developed so far. One can see that the Riemannian metric gives much less weight to large tensors, thus providing a regularization which is more robust to outliers. The anisotropic filtering further improves the results by preserving the discontinuities of the tensor scale (e.g. at the boundary of the ventricles), but also the discontinuities of the tensor orientation, which is exactly what is needed for fiber tracking in DTI.

Regularization and restoration of tensor fields

The pure diffusion is efficient to reduce the noise in the data, but it also reduces the amount of information. Moreover, the amount of smoothing is controlled by the time of diffusion (time step ε times the number of iterations), which is not an easy parameter to tune. At an infinite diffusion time, the tensor field will be completely homogeneous (or homogeneous by part for some anisotropic diffusion schemes), with a value corresponding to the mean of the measurements over the region (with Neumann boundary conditions). Thus, the absolute minimum of our regularization criterion alone is of little interest.

To keep close to the measured tensor field Σ 0 (x) while still regularizing, a more theoretically grounded approach is to consider an optimization problem with a competition between a data attachment term and a possibly non-linear anisotropic regularization term:

C(Σ) = Sim(Σ, Σ 0 ) + λ Reg(Σ).
Like before, the intrinsic evolution equation leading to a local minimum is:

Σ t+1 (x) = exp Σt(x) (-ε (∇Sim(Σ, Σ 0 ) + λ ∇Reg(Σ)(x))) .

The regularization term

As we saw in the previous section, the simplest regularization criterion is the norm of the gradient of the field Reg(F ) = Ω ∇F (x) 2 dx. To preserve the discontinuities, the gradient of this criterion (the Laplacian) may be tailored to prevent the smoothing across them, as we have done in Section 5.4. However, there is no more convergence guarantee, since this anisotropic regularization "force" may not derive from a well-posed criterion (energy). Following the pioneer work of [START_REF] Perona | Scale-space and edge detection using anisotropic diffusion[END_REF], there has been quite an extensive amount of work to propose well posed PDEs for the non-linear, anisotropic and non-stationary regularization of scalar and vector fields (see e.g. [Weickert, 1998, Sapiro, 2001] to cite only a few recent books). Some of these techniques were recently adapted to work on matrix valued fields [START_REF] Weickert | Diffusion and regularization of vector-and matrix-valued images[END_REF] (with the flat metric) or on the rotation manifolds [START_REF] Tschumperlé | Orthonormal vector sets regularization with PDE's and applications[END_REF].

One of the main idea is to replace the usual simple regularization term Reg(F ) =

Ω ∇F (x) 2 dx by an increasing function Φ of the norm of the spatial gradient: Reg(F ) = Ω Φ( ∇F (x) )dx. With some regularity conditions on the Φ-function [Aubert and Kornprobst, 2001], one can recompute the previous derivations with this Φ-function, and we end-up with:

∇Reg(F )(x) = -div Φ ( ∇F ) ∇F ∇F = - d i=1 ∂ i Φ ( ∇F ) ∇F ∂ i F .
This scheme was used in [START_REF] Chefd | Regularizing flows for constrained matrix-valued images[END_REF] with the flat Euclidean metric on tensors, in conjunction with a geometric numerical integration scheme that preserves the rank of the matrix. Their conclusion was that the rank/signature preserving flow tends to blend the orientation and diffusivity features (eigenvalue swelling effect). This ranksignature preserving flow is based on the matrix exponential but does not make any reference to a specific metric. Reformulated in our notations, their evolution equation is: Σ(x, t + dt) = exp(dt A(x, t)) Σ(x, t) exp(dt A(x, t)), where A(x, t) is implicitly related to the driving tangent vector field using: ∂ t Σ(x, t) = -W (x, t) = A(x, t)Σ(x, t) + Σ(x, t)A(x, t). From the first expression of our geodesics (Eq. 3 in Section 3.3), we can see that this is a geodesic marching scheme for our metric. However, they use the flat Euclidean metric on coefficients to evaluate the norm of the spatial gradient of the tensor field. We claim that a consistent framework should use the invariant metric. This leads to very different driving gradient fields W (x, t). Indeed large eigenvalues naturally dominate small ones with the flat metric on coefficients, which means that small eigenvalues are completely smoothed out. This produces the eigenvalue swelling effect they observed. On the contrary, the use of the invariant metric perfectly respects the discontinuities of small and large eigenvalues in our experiments of Sections 5 (even if the anisotropic diffusion PDE is not exactly the same).

We are currently investigating how to adapt the Φ-function formalism to our Riemannian tensor framework. The gradient of the modified criterion can be computed with the invariant metric like in Appendix A.1 to obtain a weighted manifold Laplacian with an additional anisotropic correction term. However, designing an efficient discrete computation scheme is more difficult. We may compute the directional derivatives using finite differences in the flat matrix space and use the intrinsic evolution scheme, but we believe that there are more efficient ways to do it using the exponential map. In the following, we keep the isotropic regularization based on the squared amplitude of the gradient.

A least-squares attachment term

Usually, one considers that the data (e.g. a scalar image or a displacement vector field F 0 (x)) are corrupted by a uniform (isotropic) Gaussian noise independent at each spatial position. With a maximum likelihood approach, this amounts to considering a least-squares criterion Sim(F ) = Ω F (x) -F 0 (x) 2 dx. Like in the previous section, we compute the first order variation by writing the Taylor expansion

Sim(F + ε H) = Sim(F ) + 2 ε Ω H(x) | F (x) -F 0 (x) dx + O(ε 2 ).
This time, the directional derivative ∂ H Sim(F ) is directly expressed using a scalar product with H in the proper functional space, so that the steepest ascent direction is ∇Sim(F ) = 2 (F (x) -F 0 (x)).

On the tensor manifold, assuming a uniform (generalized) Gaussian noise independent at each position also leads to a least-squares criterion through a maximum likelihood approach. The only difference is that it uses our Riemannian distance:

Sim(Σ) = Ω dist 2 (Σ(x) , Σ 0 (x)) dx = Ω -------→ Σ(x)Σ 0 (x) 2 Σ(x) dx.
Thanks to the properties of the exponential map, one can show that the gradient of the squared distance is: Pennec, 2004]. One can verify that this is a tangent vector at Σ whereas --→ Σ 0 Σ is not. Finally, we obtain a steepest ascent direction of our criterion which is very close to the vector case:

∇ Σ dist 2 (Σ , Σ 0 ) = -2 --→ ΣΣ 0 [
∇Sim(Σ)(x) = -2 -------→ Σ(x)Σ 0 (x). ( 15
)
6.3 A least-squares attachment term for sparsely distributed tensors

Now, let us consider the case where we do not have a dense measure of our tensor field, but only N measures Σ i at irregularly distributed sample points x i . Assuming a uniform Gaussian noise independent at each position still leads to a least-squares criterion:

Sim(Σ) = N i=1 dist 2 (Σ(x i ) , Σ i ) = Ω N i=1 dist 2 (Σ(x) , Σ i ) δ(x -x i ) dx.
In this criterion, the tensor field Σ(x) is related to the data only at the measurement points x i through the Dirac distributions δ(x -x i ). If the introduction of distributions may be dealt with for the theoretical differentiation of the criterion with respect to the continuous tensor field Σ, it is a real problem for the numerical implementation. In order to regularize the problem, we consider the Dirac distribution as the limit of the Gaussian function G σ when σ goes to zero. Using that scheme, our criterion becomes the limit case σ = 0 of:

Sim σ (Σ) = Ω N i=1 dist 2 (Σ(x) , Σ i ) G σ (x -x i ) dx. (16)
From a practical point of view, we need to use a value of σ which is of the order of the spatial resolution of the grid on which Σ(x) is evaluated, so that all measures can at least influence the neighboring nodes. Now that we came back to a smooth criterion, we may differentiate it exactly as we did for the dense measurement setup. The first order variation is:

Sim σ (Σ + εW ) = Sim σ (Σ) -2 ε N i=1 Ω W (x) G σ (x -x i ) ----→ Σ(x)Σ i dx + O(ε 2 ),
so that we get:

∇Sim σ (x) = -2 N i=1 G σ (x -x i ) ----→ Σ(x)Σ i . (17)

Interpolation through diffusion

With the sparse data attachment term ( 16) and the isotropic first order regularization term (10), we are looking for a tensor field that minimizes its spatial variations while interpolating (or more precisely approximating at the desired precision) the measurement values:

C(Σ) = N i=1 Ω G σ (x -x i ) dist 2 (Σ(x i ) , Σ i ) + λ Ω ∇Σ(x) 2 Σ(x) dx.
According to the previous sections, the gradient of this criterion is

∇C(Σ)(x) = -2 N i=1 G σ (x -x i ) ----→ Σ(x)Σ i -2 λ ∆Σ(x).
Using our finite difference approximation scheme (Eq. 13), the intrinsic geodesic gradient descent scheme (Sec. 3.6) is finally:

Σ t+1 (x) = exp Σt(x) ε N i=1 G σ (x -x i ) ----→ Σ(x)Σ i + λ u∈V ---------→ Σ(x)Σ(x + u). u 2 (18)
Last but not least, we need an initialization of the tensor field Σ 0 (x) to obtain a fully operational algorithm. This is easily done with any radial basis function approximation, for instance the renormalized Gaussian scheme that we investigated in Section 4.4. Figure 10 displays the result of this algorithm on the interpolation between 4 tensors. On can see that the soft closest point approximation is well regularized into a constant field equal to the mean of the four tensors if data attachment term is neglected. On the contrary, a very small value of λ is sufficient for regularizing the field between known tensors (as soon as σ is much smaller than the typical spatial distance between two measurements).

The choice of the initialization is a critical issue from a computational point of view. For instance, starting with a constant (or any harmonic) field is a bad idea: there is a null Laplacian everywhere, except at the immediate neighborhood of the sparse tensors, exactly where the data attachment term acts. Thus, we have a potentially destructive competition between the two terms of the criterion in very localized area. On the contrary, starting with a soft closest point approximation leads to a Laplacian which is non null on the boundaries of the Voronoi cells of the measurement points, i.e. the farthest possible place from the sparse measures. In that case, the Laplacian regularization will spread from these boundaries with no constraints until it reaches the counterbalancing forces of the data attachment term in the immediate vicinity of the sparse measurements. Thus, we may expect to reach the maximal efficiency in terms of convergence rate.

Conclusion

We propose in this paper an affine invariant metric that endows the space of positive define symmetric matrices (tensors) with a very regular manifold structure. In particular, tensors with null and infinite eigenvalues are both at an infinite distance of any positive definite symmetric matrix: the cone of positive definite symmetric matrices is replaced by a space which has an infinite development in each of its n(n + 1)/2 directions. Moreover, there is one and only one geodesic joining any two tensors, and we can even define globally consistent orthonormal coordinate systems of the tangent spaces. Thus, the structure we obtain is very close to a vector space, except that the space is curved. We exemplify some the the good metric properties for some simple statistical operations. For instance, the Karcher mean in Riemannian manifolds has to be defined through a distance-based variational formulation. With our invariant metric on tensors, the existence and uniqueness is insured, which is generally not the case.

A second contribution of the paper is the application of this framework to important geometric data processing problem such as interpolation, filtering, diffusion and restoration of tensor fields. We show that interpolation and Gaussian filtering can be tackled efficiently through a weighted mean computation. However, if weights are easy to define for regularly sampled tensors (e.g. for linear to tri-linear interpolation), the problem proved to be more difficult for irregularly sampled values. The solution we propose is to consider this type of interpolation as a statistical restoration problem where we want to retrieve a regular tensor field between (possibly noisy) measured tensor values at sparse points. This type of problem is usually solved using a PDE evolution equation. We show that the usual linear regularization (minimizing the magnitude of the gradient) and some anisotropic diffusion schemes can be adapted to our Riemannian framework, provided that the metric of the tensor space is taken into account. We also provide intrinsic numerical schemes for the computation of the gradient and Laplace-Beltrami operators. Finally, simple statistical considerations led us to propose least-squares data attachment criteria for dense and sparsely distributed tensor fields. The differentiation of these criteria is particularly efficient thanks to the use of the Riemannian distance inherited from the chosen metric.

From a theoretical point of view, this paper is a striking illustration of the general framework we are developing since [Pennec, 1996] to provide a rigorous computing environment for geometric objects. This framework is based on the choice of a Riemannian metric on one side, which leads to powerful differential geometry tools such as the the exponential maps and geodesic marching techniques, and on the transformation of linear combinations or integrals into minimization problems on the other side. The Karcher mean and the generalized Gaussian distribution are a typical example that we previously investigated [Pennec, 2004]. In the present paper, we provide new examples with interpolation, filtering and PDEs on Riemannian-valued fields.

Many research avenues are still left open, in particular the choice of the metric to use. In a more practical domain, we believe that investigating new intrinsic numerical schemes to compute the derivatives in the PDEs could lead to important gains in accuracy and efficiency. Last but not least, all the results presented in this paper still need to be confronted to other existing methods and validated in the context of medical DTI applications. We are currently investigating another very interesting application field in collaboration with P. Thompson and A. Toga at UCLA: the modeling and analysis of the variability of the brain anatomy.

A Tensor regularization: the Laplace-Beltrami operator

A.1 Gradient of the L 2 regularization of a tensor field Let Σ(x) ∈ Sym + n be a tensor field over R d , and (x 1 , . . . x d ) be an orthonormal coordinate system. To simplify the notations, we use in this section ∂ i for the spatial derivative ∂/(∂x i ) and we do not specify the (spatial) integration variable x. The L 2 regularization criterion

Reg(Σ) = Ω ∇Σ(x) 2 Σ(x) dx = d i=1 Ω ∂ i Σ 2 Σ = d i=1 Ω
Tr (∂ i Σ) Σ (-1) (∂ i Σ) Σ (-1) .

Using the Taylor expansion (Σ + εW ) (-1) = Σ (-1) -εΣ (-1) W Σ (-1) + O(ε 2 ) in the the Taylor expansion of our regularization criterion and identifying the first order term to Reg(Σ + εW ) = Reg(Σ) + ε ∂ W Reg + O(ε 2 ), we get the directional derivative:

∂ W Reg = 2 d i=1 Ω Tr (∂ i W ) Σ (-1) (∂ i Σ) Σ (-1) -(∂ i Σ) Σ (-1) (∂ i Σ) Σ (-1) W Σ (-1)
The main goal is to find out the field of tangent vectors ∇Reg(x) ∈ T Σ(x) Sym + n such that, by definition of the gradient, we have the equality: ∂ W Reg = Ω W | ∇Reg Σ dx for every field of tangent vectors W (x) ∈ T Σ(x) Sym + n . As the above expression of ∂ W Reg is in the standard Euclidean chart (matrix coefficients), we shall safely use the computations of the previous sections. Notice that we are using the matrix coefficients only as a chart and not as a metric. Let Λ i = Σ (-1) (∂ i Σ) Σ (-1) . We get:

∂ W Reg = 2 d i=1 Ω Tr((∂ i W ) Λ i ) -W | (∂ i Σ) Σ (-1) (∂ i Σ) Σ
Now, assuming the proper Neumann boundary conditions, we can apply the previous integration by part formula Ω Tr((∂ i W ) Λ i ) = -Ω Tr(W (∂ i Λ i )) to the first term:

∂ W Reg = -2 d i=1 Ω Tr W Σ (-1) (Σ(∂ i Λ i )Σ) Σ (-1) + W | (∂ i Σ) Σ (-1) (∂ i Σ) Σ = -2 d i=1 Ω W | Σ(∂ i Λ i )Σ + (∂ i Σ) Σ (-1) (∂ i Σ) Σ
We have obtained the expression that defines the gradient of our regularization criterion:

∇Reg = -2 d i=1 Σ (∂ i Λ i ) Σ + (∂ i Σ) Σ (-1) (∂ i Σ)
To compute explicitly its value, let us observe first that ∂ i Σ (-1) = -Σ (-1) (∂ i Σ) Σ (-1) because ∂ i (Σ (-1) Σ) = 0. Thus, thanks to the chain rule, we have:

Σ (∂ i Λ i ) Σ = Σ ∂ i (Σ (-1) (∂ i Σ) Σ (-1) ) Σ = ∂ 2 i Σ -2 (∂ i Σ) Σ (-1) (∂ i Σ)

Figure 2 :

 2 Figure 2: Mean of random Gaussian tensors. Left: Typical evolution of the distance between successive iterations of the mean computation. The convergence is clearly very fast. Right: Histogram of the renormalized Mahalanobis distance µ 2 = N dist( Σ, Σ) 2 /γ 2 between the computed and the exact mean tensors. The curve is the pdf of the χ 2 6 distribution.

Figure 3 :

 3 Figure 3: Top: Linear interpolation between 2D tensors of eigenvalues (5,1) horizontally and (1,50) at 45 degrees. Left: interpolation in the standard matrix space (interpolation of the coefficients). Right: geodesic interpolation in our Riemannian space. Bottom: evolution of the eigenvalues, their mean (i.e. trace of the matrix) and product (i.e. determinant of the matrix or volume of the ellipsoid).

Figure 4 :

 4 Figure 4: Top left: Bi-linear interpolation between the four 2D tensors at the corners in the standard matrix space (interpolation of the coefficients). Top right: Equivalent bi-linear interpolation in our Riemannian space. Bottom left: A slice of the tri-linear interpolation between 3D tensors in the standard matrix space (interpolation of the coefficients). Bottom right: Equivalent tri-linear interpolation in our Riemannian space.

Figure 5 :

 5 Figure 5: Regularization of a DTI slice around the corpus callosum by isotropic Gaussian filtering.

Figure 6 :

 6 Figure 6: Norm of the gradient of the tensor field. Left: computed on the coefficients with Eq. 7

Figure 7 :

 7 Figure 7: Left: 3D synthetic tensor field with a clear discontinuity. Middle: The field has been corrupted by a Gaussian noise (in the Riemannian sense). Right: result of the regularization after 30 iterations (time step ε = 0.01).

Figure 8 :

 8 Figure 8: Anisotropic filtering of a DTI slice (time step 0.01, κ = 0.046). From left to right: at the beginning, after 10 and after 50 iterations. Top: A 3D view of the tensors as ellipsoids. The color codes for the direction of the principal eigenvector. The results could be compared with the isotropic Gaussian filtering displayed in Figure 5. Middle: Riemannian norm of the gradient. Bottom: fractional anisotropy.

Figure 9 :

 9 Figure 9: Closeup on the results of the different filtering methods around the splenium of the corpus callosum. The color codes for the direction of the principal eigenvector (red: left-right, green: posterior-anterior, blue: inferior-superior). Top left: Original image. Top right: Gaussian filtering using the flat metric (5x5 window, σ = 2.0). This metric gives too much weight to tensors with large eigenvalues, thus leading to clear outliers in the ventricles or in the middle of the splenium tract. Bottom left: Gaussian filtering using the Riemannian metric (5x5 window, σ = 2.0). Outliers disappeared, but the discontinuities are not well preserved, for instance in the ventricles at the level of the cortico-spinal tracts (upper-middle part of the images). Bottom right: Anisotropic filtering in the Riemannian framework (time step 0.01, 50 iterations). The ventricles boundary is very well conserved with an anisotropic filter and both isotropic (ventricles) and anisotropic (splenium) regions are regularized. Note that the U-shaped tracts at the boundary of the grey/white matter (lower left and right corners of each image) are preserved with an anisotropic filter and not with a Gaussian filter.

Figure 10 :

 10 Figure 10: Interpolation and extrapolation of tensor values from four measurements using diffusion. Top left: The four initial tensor measurements. Top right: Initialization of the tensor field using a soft closest point interpolation (mean of the four tensors with a renormalized spatial Gaussian influence). Bottom left: result of the diffusion without the data attachment term (1000 iterations, time-step ε = 1, λ = +∞). Bottom right: result of the diffusion with an attachment term after (1000 iterations, time-step ε = 1, λ = 0.01, σ = 1 pixel of the reconstruction grid). The algorithm did in fact converge in about 100 iterations.

To be mathematically correct, we should consider the quotient space Sym + n = GL + n /SO n instead of Sym + n = GL n /O n so that all spaces are simply connected.

On most homogeneous manifolds, this can only be realized locally. For instance, on the sphere, there is a singularity at the antipodal point of the chosen origin for any otherwise smooth placement function.

Eventually, we end up with ∇Reg(Σ) = -2∆Σ, where ∆ is the Laplace-Beltrami operator on our manifold:

A.2 Numerical implementation of the Laplace-Beltrami operator

From the Taylor expansion of a tensor field (considered as a matrix field) Σ at x, we have

). Thus, we may approximate the first and second order tensor derivatives by their Euclidean derivatives:

This finally gives us a fourth order approximation of the Laplace-Beltrami operator in the spatial direction u:

However, this numerical scheme is extrinsic since it is based on (Euclidean) differences of tensors. We propose here an intrinsic scheme based on the exponential chart at the current point: we claim that ---------→ Σ(x)Σ(x + u) + ---------→ Σ(x)Σ(x -u) is a forth order approximation of the Laplace Beltrami operator in the direction u. Indeed, we have

where we put

From the Log series (Eq. 1), we get:

The Taylor expansion of ----------→ Σ(x)Σ(x -εu) is obtained by replacing ε by -ε, so that we finally end up with -