Emmanuel Chevallier

Emmanuel Chevallier Morphology

Directeur De

Jesús Angulo

Jury M Fernand

Professeur M Michel Berthier

Professeur Mme Isabelle Bloch

Professeur Rapporteur Mme Gabriele Steidl

Senior Scientist, THALES Rapporteur M Frédéric Barbaresco

Examinateur M Jasper Van De Gronde

Professeur M Christian Ronse

Morphology, Geometry and Statistics in non-standard imaging

Keywords: Morphologie Mathématique, Ordres Totaux, Traitement d'Image, Estimation de Densités, Géometrie Riemannienne, Métrique de Fisher, Métrique de Wasserstein, Histogrammes d Mathematical Morphology, Total Orders, Image Processing, Density Estimation, Riemanninan Geometry, Fisher Metric, Wasserstein Metric, Image Histograms iii Mathematical Morphology, Total Orders, Image Processing, Density Estimation, Riemanninan Geometry, Fisher Metric, Wasserstein Metric, Image Histograms

Le traitement d'images numériques a suivi l'évolution de l'électronique et de l'informatique.

Il est maintenant courant de manipuler des images à valeur non pas dans {0, 1}, mais dans des variétés ou des distributions de probabilités. C'est le cas par exemple des images couleurs où de l'imagerie du tenseur de diusion (DTI). Chaque type d'image possède ses propres structures algèbriques, topologiques et géometriques. Ainsi, les techniques existantes de traitement d'image doivent être adaptés lorsqu'elles sont appliquées à de nouvelles modalités d'imagerie. Lorsque l'on maniuple de nouveaux types d'espaces de valeurs, les précédant opérateurs peuvent rarement être utilisés tel quel. Même si les notions sous-jacentes ont encore un sens, un travail doit être mené an de les exprimer dans le nouveau contexte. Cette thèse est composée de deux parties indépendantes. La première, Morphologie matématiques pour les images non standards , concerne l'extension de la morphologie mathématique à des cas particuliers où l'espace des valeurs de l'image ne possède pas de structure d'ordre canonique. Le chapitre 2 formalise et démontre le problème de l'irregularité des ordres totaux dans les espaces métriques. Le résultat principal de ce chapitre montre qu'étant donné un ordre total dans un espace vectoriel multidimensionnel, il existe toujours des images à valeur dans cet espace tel que les dilatations et les érosions morphologiques soient irrégulières et incohérentes. Le chapitre 3 est une tentative d'extension de la morphologie mathématique aux images à valeur dans un ensemble de labels non ordonnés. La deuxième partie de la thèse, Estimation de densités de probabilités dans les espaces de Riemann concerne l'adaptation des techniques classiques d'estimation de densités non paramétriques à certaines variétés Riemanniennes. Le chapitre 5 est un travail sur les histogrammes d'images couleur dans le cadre de métriques perceptuelles. L'idée principale de ce chapitre consiste à calculer les histogrammes suivant une approximation euclidienne local de la métrique perceptuelle, et non une approximation globale comme dans les espaces perceptuels standards. Le chapitre 6 est une étude sur l'estimation de densité lorsque les données sont des lois Gaussiennes. Diérentes techniques y sont analysées. Le résultat principal est l'expression de noyaux pour la métrique de Wasserstein.

Introduction Résumé

Ce chapitre présente les motivations et les thèmes abordés dans cette thèse. L'augmentation progressive de la mémoire disponible et l'évolution des dispositifs physiques d'acquisition ont permis l'apparition d'images de plus en plus complexes. Les techniques existantes de traitement doivent être adaptées à chaque nouvelle modalité d'image. Deux problématiques sont abordées dans cette thèse. La première concerne la morphologie mathématique, dans des cas particuliers où l'espace des valeurs de l'image ne dispose pas intrinsèquement d'une structure de treillis. La seconde concerne l'estimation de densité dans des variétés Riemanniennes. Ce problème apparaît notamment lors de la construction d'histogrammes d'images à valeurs dans des variétés Riemanniennes.

Motivation: Non standard images

Since its origins, digital signal processing follows the evolution of electronic technology. In the sixties, with the development of computers, image processing has become a discipline in itself. For technological reasons, rst computers could only handle images composed of two distinct levels, 0 and 1. Thus, the rst techniques of image processing including mathematical morphology where developed for binary images. Since the sixties, the development of sensors and computers has widened digital image processing to numerous new forms of images. The nature of images evolves on two levels. On the one hand, the spatial support has evolved from a regular grid of pixels to an arbitrary graph. On the other hand, the value contained in each pixel or node of the graph takes more and more complex forms.

The increase of available memory has rst enabled the transition from binary to gray-scale images. The memory is now sucient to describe a signicant amount of manifolds such as colors or probability distributions. Diusion tensor imaging (DTI) is a typical case of non standard images. A physical device measures the motion of water molecules and models it by a centered 3-dimensional Gaussian distribution. The set of 3-dimensional Gaussian distributions are a 6-dimensional manifold, naturally represented by the 3 × 3 symmetric positive denite matrices. Another case of images valued in distributions arise from the popularization of superpixels. In order to simplify images, neighbor pixels can be gathered according to a given criterion, to form a set of superpixels. Usually a single value is chosen to represent a superpixel. However, it is sometimes interesting to preserve more information and to assign a distribution of values to each superpixel instead of a single value.

Each kind of images has its own algebraic, topological and geometric properties. These properties depend on the graph of pixels, on the value space, and on the nature of the observed objects. In this thesis, we leave aside the problem of the graph of pixels and focus on the value space and its structure. Usually, an image processing operator depends on the structure of the image. Some of them have a certain generality. For instance, the 1. Introduction mean and the median lters are dened on a large variety of spaces. However the denition of the lters depends on the value space of the image: the mean and the median, originally dened using addition, division and an order, have been reformulated in metric spaces.

When dealing with new kind of value spaces, former operators can rarely be used as they are. Even if the underlying notion has still a meaning, a work must be carried out in order to express it in the new context.

Statistics in non-Euclidean spaces

The adaptation of standard quantities of statistic to new spaces is an important problem.

For instance, the expected value was rst dened using the algebraic structure of R. This denition is naturally extended to random phenomena valued in R n . Consider a distribution µ, taking values on the surface of the earth S, for instance the distribution of the population in a country. The notion of mean or expected value has an intuitive meaning but can not be dened from algebraic operations. One has to nd the counterpart of the previous notion based on the only intrinsic structure of S: the notion of distance d(x ∈ S, x ∈ S). Without entering uniqueness problems, the following quantities

E S = argmin x x ∈S d(x, x ) 2 dµ, E S emp (n) = argmin x i d(x, x i ) 2 .
correspond to the original expected value and empirical expected value in R when replacing the distance on S by the distance on R. Furthermore, they verify similar convergence

properties (E S emp (n) → E S ) for many standard distances d. A mean formalized according to this denition is called a Fréchet mean. The situation is similar for the space of colors.

The set of colors has intrinsically no algebraic structure, but posses intrinsic metrics. Thus a mean between colors should be a Fréchet mean.

In this thesis we are not specically interested in the notion of mean but in the notion of density of a distribution. Unlike the mean where the denition must be adapted, the denition of the density is not a problem as long as the value space has a reference measure.

However the estimation of a probability density from a set of samples strongly depends on the context. Consider a space X and two metrics d 1 and d 2 on X . The most straightforward way of estimating a density is to use a tilling of the space. Estimating a density consists in counting a number of samples in each bin of the tilling and dividing by their surface.

Finding good tilling is a dicult problem that has to be studied for each case (X , d i ).

The same holds for every density estimation techniques. Chapter 6 studies deeply density estimation on spaces of Gaussian and multivariate Gaussian distributions.

Mathematical morphology for non standard images

Similarly to statistics, it is possible to observe how mathematical morphology have progressed with the evolution of the value spaces of images. Mathematical morphology is one of the rst sound theories of digital image processing. When dealing with binary images, three neighborhood operators naturally arise: replacing the value of a pixel by the majority, the minimal or the maximal value in the neighborhood. Replacing the value of a pixel by the majority value in a neighborhood is called a median lter. The two other operators called erosion, ε, and dilation, δ, respectively form the famous couple of operators at the origin of mathematical morphology. This pair of operators present remarkable stability properties. For a suitable choice of neighborhood, the operator formed by the composition of an erosion and a dilation is stable after one iteration, that is to say ε

• δ = ε • δ • ε • δ.
By symmetry, or by duality in the morphological vocabulary, the same holds for δ • ε. The composition ε • δ is called a closing and δ • ε an opening. This stability property, also called idempotency, together with the scale-space semigroups of openings and closings, are at the heart of morphological image ltering. The erosion and the dilation can be extended to grey-scale images by taking the minimal and maximal value in the neighborhood according to the standard order on reals. This generalization preserves the fundamental properties of stability. The theory is now generalized in terms of complete lattice theory, a lattice being a partially ordered set where couples of elements have a unique supremum and inmum.

Pairs of erosion and dilation verifying a pseudo-inverse property called adjunction generate idempotent compositions. In this framework, the value space of the image has to be endowed with a complete lattice structure in order to use morphological operators. This is natural for grey-scale images but often problematic for images valued in multidimensional spaces. Colors can be represented in a three dimensional vector space. Standard lattice structures used to process colors all have important drawbacks. Partial orders tend to introduce false colors while total orders are not compatible with the Euclidean distance as

Chapter 2 shows.

The idea of formulating the best total order for a given image is the object of Chapter 2.

Several parallel works to this thesis aim at dening morphological operators without the notion of complete lattice. Chapter 3 proposes a morphological framework not based on lattices, for images valued in a set of independent categories. Approach introduced by [START_REF] Gronde | Sponges for Generalized Morphology[END_REF] proposes to drop the associative law of the lattice structure.

Unlike the proposition of chapter 3 their framework enables to work in a continuous space, however the idempotence of openings and closings no longer holds. Meanwhile, work by [START_REF] Cannon | A Color Tree of Shapes with Illustrations on Filtering, Simplication, and Segmentation[END_REF] build a tree of shapes on multivariate images without requiring an order structure.

Thesis organization and main results

This thesis was originally focused on mathematical morphology. In the chapter Unsupervised Morphology" of the PhD thesis [START_REF] Velasco-Forero | Topics in mathematical morphology for multivariate images[END_REF], the author proposes to build an order on the value space of an image from the distribution of the values using a notion of statistical depth. In theory, there is no need to compute a density associated to the point cloud of values in order to compute the statistical depth. Indeed, the statistical depth is dened directly from the point cloud of values. However estimating the density of the distribution of the values has an interest in term of computation time of the statistical depth. Trying to apply this idea to order DT I data, we encountered the problem of density estimation of symmetric positive denite matrices. This problem took a growing importance and gradually led to an autonomous part of the thesis on density estimation on Riemannian spaces. Meanwhile our development of mathematical morphology changed direction and the original idea of ordering DT I data according to a statistical depth has been abandoned.

The thesis is thus composed of two independent parts. The rst part named Mathe-matical morphology on non-standard images is composed of two chapters:

Chapter 2 formalizes and demonstrates the irregularity issue of total orders in metric spaces, its main result being the lemma 2.3.1. This problem has been known for a while in the mathematical morphology community, see pioneering work by [START_REF] Chanussot | Total Ordering Based on Space Filling Curves for Multivalued Morphology[END_REF] or more recently, by [START_REF] Flórez-Revuelta | Ordering of the RGB Space with a Growing Self-organizing Network. Application to Color Mathematical Morphology[END_REF]. The problem is reformulated so as to better take into account the needs of mathematical morphology.

Our original formal result states that any total order on a multidimensional space possesses a form of irregularity that generates inconsistent and irregular erosions and dilations. A method is proposed to build total orders adapted to each image, by minimizing the irregularity of the order.

n-ary images refer in the thesis to images valued in a set of n elements or categories, with no specic mathematical structure. Each value is called a label. n-ary images appear when the studied objects have equivalent semantic levels. Consider for instance a trac image with blue, red and yellow cars. Blue, red and yellow represent objects of the same semantic level. In this situation, distances and orders on colors present an arbitrary character: saying that blue is closer to yellow than red or saying that blue<yellow<red has no meaning in terms of the studied objects. This situation often appears in pixelwise classied images. It is dicult to address data that are intrinsically unordered using standard morphological operators. Chapter 3 is an attempt to generalize morphology on binary images to discrete n-ary images.

The previous considerations already have motivated several works, see [START_REF] Busch | Morphological Operations for Color-Coded Images[END_REF] and [START_REF] Ronse | Morphology on Label Images: Flat-Type Operators and Connections[END_REF], but led to dierent directions. The main originality of this work is the introduction of multiple erosions and dilations. The main results of the chapter are,

The denition of an erosion ε i and a dilation δ i per label; A proof of the idempotence of γ i = ε i • δ i and ϕ i = δ i • ε i . Recall that this results is standard in the binary case by adjunction theory;

A proof of semigroup ε i,R 1 • ε i,R 2 = ε i,R 1 +R 2 .
In other words, an erosion with a ball of radius R 1 composed by a erosion with a ball of radius R 2 is equivalent to an erosion with a ball of radius R 1 + R 2 . The result also holds for dilation but is directly deducted from the binary case. This enables to express the n-ary operations as ow propagations or operator composition;

The operator F = (ϕ 1 • ... • ϕ n ), involving dierent labels, is no longer idempotent. However, F k = (ϕ 1 • ... • ϕ n ) k is stable when k goes to innity. F is a good candidate to lter out image objects with shape/shape criteria;

A rst extension to the continuous n-ary case.

The second part entitled Probability density estimation on Riemannian spaces is composed of three chapters:

Chapter 4 introduces a background on probability and geometry for a typical reader from image processing community.

Chapter 5 addresses the issue of image histogram computation. A formalization of the notion of image histogram in term of density estimation is rst proposed. The second part of the chapter is a study on perceptual color image histograms. The perceptual metric on colors is assumed to be Riemannian. Local scalar products are provided by standard ellipse or ellipsoid datasets, such as the MacAdam ellipses. Usually perceptual histograms are built according to the Euclidean structure of perceptual spaces such as the Lab space. The use of local Euclidean approximations enables to follow more precisely the Riemannian metric without increasing the computational complexity.

Chapter 6 is a guide to non parametric probability density estimation on spaces of multivariate Gaussian distributions. In other words it provides methods to estimate a density when the set of samples {x i } is a set of Gaussian distributions. This situation is more and more common in signal and image processing. The case of Gaussian distributions also includes the space of positive denite matrices and the hyperbolic space of dimension two. The construction of the estimated density depends on the choice of a metric. The metrics studied in this chapter are the main standard metrics on Gaussian distributions: the Fisher metric, the Wasserstein metric and Euclidean metrics on parameters. The main results are, An expression of kernels for kernel density estimation on multivariate Gaussian distributions when the space is endowed with the Wasserstein metric;

An expression of kernels for kernel density estimation on multivariate centered Gaussian distributions when the space is endowed with the Fisher metric. This result is not entirely original, the dierent factors composing the kernel are already present in the specialized literature, but are here systematically developed;

The study of partial quantities (mean, eigenvalues and rotation) under the different metrics.

Part I.

Mathematical morphology on non standard images 2. The irregularity issue of total orders on metric spaces and its consequences for mathematical morphology

Résumé

Nous nous intéressons dans ce chapitre au problème de l'irrégularité des ordres totaux dans les espaces métriques pour la morphologie mathématique. Nous donnons tout d'abord une formulation rigoureuse du problème. Une nouvelle approche est proposée pour y faire face, en adaptant l'ordre à l'image à traiter. Étant donné une image et un ordre total, nous dénissons une fonction de coût évaluant l'importance des conits lors des opérateurs morphologiques. L'ordre proposé est alors une minimisation de cette fonction de coût. Un des atouts de cette méthode est sa généralité : le seul ingrédient nécessaire à la construction de l'ordre est un graphe des distances entre les valeurs présentes dans l'image. L'ordre peut donc être calculé pour toute les images à valeur dans des espaces métriques où l'expression des distances est connue. Nous présentons quelques résultats sur des images couleurs, sur l'imagerie du tenseur de diusion (DTI), et sur des images à valeur dans le demi plan supérieur de Poincaré.

Ce chapitre est principalement adapté de Chevallier and Angulo (2015).

Abstract

We address in this chapter the problem of irregularity of total orders in metric spaces and its implications for mathematical morphology. We rst give a rigorous formulation of the problem. Then, a new approach is proposed to tackle the issue by adapting the order to the image to be processed. Given an image and a total order, we dene a cost that evaluates the importance of the conict for morphological processing. The proposed order is then built as a minimization of this cost function. One of the strength of the proposed framework is its generality: the only ingredient required to build the total order is the graph of distances between values of the image. The adapted order can be computed for any image valued in a metric space where the distance is explicitly known. We present results for color images, diusion tensor images (DTI) and images valued in the hyperbolic upper half-plane.

This chapter is mainly adapted from Chevallier and Angulo (2015).

Introduction

Since its emergence in the sixties, mathematical morphology has become one of the major theory of nonlinear image processing. Originally used for binary images in [START_REF] Matheron | Random sets and integral geometry[END_REF], the theory has followed the technical evolution of computer science which has enabled the manipulation of more and more complex images, see [START_REF] Serra | Image analysis and mathematical morphology[END_REF]. The set theory was sucient to study binary images. Later, the emergence of gray-scale images required the introduction of the notion of order. From the works by [START_REF] Serra | Image Analysis and Mathematical Morphology[END_REF] and [START_REF] Heijmans | The algebraic basis of mathematical morphology I. Dilations and erosions[END_REF], the theory of mathematical morphology is now fully based on the lattice theory. Value spaces endowed with a total order form the most comfortable framework for morphological processing. However, if it is natural to endow gray-scale images with a total order, the task is more dicult when the pixel values do not have a one-dimensional structure, see [START_REF] Serra | Anamorphoses and function lattices[END_REF] for a pioneering discussion on the topic. Indeed, we show that the information contained in a total order is too weak to completely represent the value space. In many situations, the use of a partial order such as a product order on a vector space is preferred. Using a product order is equivalent to processing components independently. Then this order structure becomes natural but some information about the geometry of the original value space is lost. Both choices present a loss of information.

The interest of product orderings in mathematical morphology is still a matter of recent research, see for instance [START_REF] Gronde | Group-Invariant Frames for Colour Morphology[END_REF] and [START_REF] Burgeth | Morphology for Color Images via Loewner Order for Matrix Fields[END_REF], mainly by considering the geometric and invariant properties of the underlying space.

Other recent works on partial ordering for morphological operators on vector images were motivated by the need of taking into account some prior information to order vectors:

either by learning the order from training samples in Velasco-Forero and [START_REF] Velasco-Forero | Supervised Ordering in : Application to Morphological Processing of Hyperspectral Images[END_REF] or by building the order according to the outlier distribution in Velasco-Forero and Angulo (2012). Therefore, partial ordering can become image adaptive and consequently leads to more relevant morphological operators.

On the other hand, the study of total orders remains mainly limited to lexicographic orders, see [START_REF] Hanbury | Morphological operators on the unit circle[END_REF]; [START_REF] Angulo | Morphological colour operators in totally ordered lattices based on distances: Application to image ltering, enhancement and analysis[END_REF][START_REF] Angulo | Geometric algebra colour image representations and derived total orderings for morphological operators Part I: Colour quaternions[END_REF]; [START_REF] Aptoula | A Comparative Study on Multivariate Mathematical Morphology[END_REF].

Since the emergence of multivariate images, very few papers have addressed the problem of total orders in a general way. Previous approaches of total ordering mainly focused on building regular orders on hypercubes see [START_REF] Chanussot | Total Ordering Based on Space Filling Curves for Multivalued Morphology[END_REF] and [START_REF] Flórez-Revuelta | Ordering of the RGB Space with a Growing Self-organizing Network. Application to Color Mathematical Morphology[END_REF]. Beyond the theoretical result, the two ideas of the present work are the following. Firstly, given an image we restrict the value space to the values that are actually present in the image. This idea is already present in [START_REF] Lezoray | Graph-based ordering scheme for color image ltering[END_REF] where they order the values present in the image according to a Hamiltonian path. Although the motivation is dierent, the idea is also present in [START_REF] Ledda | Majority Ordering for Colour Mathematical Morphology[END_REF]. Secondly, we take into account the locations of the values in the image. Given an image, these considerations enable to nd more regular orders. One can therefore consider our approach as an image adapted total ordering. We formulate this task as an optimization problem which cannot be solved using classical optimization techniques. Thus, we introduce a hierarchical recursive algorithm aiming at nding an approximated solution.

The chapter is organized as follows.

Section 2.2 provides a background on notations and basic notions.

A discussion on existing ordering strategies on vector spaces is given in Section 2.3.

In particular, a strong theoretical result points out the limitation of total orderings in terms of discontinuity of vector morphological operators.

Section 2.4 gives the motivation for introducing an order adapted to a given image.

A cost function is introduced to measure the quality of a total order regarding the topological conict.

Notations and recalls

The recursive algorithm developed to minimize the cost function is fully described in Section 2.5.

A short discussion on the invariance properties of the adapted order is given in Section 2.6.

Section 2.7 presents some results of morphological image processing using the image adapted total order. We present results for color images, diusion tensor images (DTI) and images valued in the hyperbolic upper half-plane.

Conclusions and perspectives close the chapter in Section 2.8.

Notations and recalls

We set here a few notations and recall elementary operators of mathematical morphology, see for instance [START_REF] Serra | Image analysis and mathematical morphology[END_REF]; [START_REF] Soille | Morphological Image Analysis[END_REF]. Let us consider an image I as a function:

I : Ω → V p → I(p)
where Ω is the support space of pixels p: typically Ω ⊂ Z 2 or Z 3 for discrete images. The pixel values of the image belong to the space V. Typically we have V ⊂ R for grey-scale images, V ⊂ R n for multivariate vector images, or V ⊂ M for manifold valued images.

We address here images where V is any metric space. Points in V will generally be called colors. We denote by I(Ω) ⊂ V the set of colors of V present in the image I(p).

Unlike linear processing mainly based on linear convolution (i.e., weighted average), mathematical morphology is based on sup and inf-convolution. The numerous choices of convolution kernels lead to a large variety of processing. Thus, the two basic operators of mathematical morphology are the erosion and the dilation of an image I(p), I : Ω → R, by B ⊂ Ω given respectively by

ε B (I)(p) = inf q∈B(p) {I(q)} , (2.1) δ B (I)(p) = sup q∈ B(p) {I(q)} , (2.2)
where the set B denes the structuring element (the equivalent of the convolution kernel), B is the transpose of B (i.e., symmetric set with respect to the origin), and B(p) denes the neighborhood of p according to the shape of B. Note that here we only focus on at structuring elements. Other morphological lters, such as the opening γ B (I) and closing ϕ B (I), are obtained by composition of dilation and erosion; i.e,

γ B (I) = δ B (ε B (I)) , (2.3) ϕ B (I) = ε B (δ B (I)) .
(2.4)

More evolved lters and transforms are obtained from composition of openings/closings.

Another nonlinear operator particularly useful in image denoising, also based on ordering, is the median lter:

m B (I)(p) = med q∈B(p)
{I(q)} .

(2.5) Figure 2.1.: The metric space X is here a rectangle of the Euclidean plan

Existing total orders

The problem of total ordering for multivariate images is a relatively well known problem in mathematical morphology. The essential diculty is that the topology induced by a total order on a multidimensional space cannot reproduce the natural topology of the vector space, see discussion in [START_REF] Chanussot | Approches vectorielles ou marginales pour le traitement d'images multi-composantes[END_REF]. Arising as a milestone limitation, we have the following theoretical result.

Lemma 2.3.1 Let (X, d) be a metric space endowed with a total order ≤. Suppose that there exist a positive real number R and three points x 1 , x 2 , x 3 ∈ X such that d(x i , x j ) > R for i = j and that the complementary B C (x i , R) of each ball B(x i , R) is connected, as in Fig. 2.1. Then for all r > 0, there exist three points a, b and c in X such that

   a ≤ b ≤ c, d(a, b) ≥ R, d(a, c) ≤ r.
Proof We can assume that x 1 < x 2 < x 3 . We argue by contradiction and assume that there exists r > 0 such that for all a, b, c ∈ X one at least of the three above conditions doesn't hold. The idea of the proof is to show that under this assumption, we have either

∀y / ∈ B(x 2 , R), y < x 2 or ∀y / ∈ B(x 2 , R), y > x 2 .
This would contradict the existence of the triple (x 1 , x 2 , x 3 ) with x 1 < x 2 < x 3 . By our assumption, for any triple (a, b, c) of points in X, the three inequalities a ≤ b, d(a, b) ≥ R and d(a, c) ≤ r imply c < b. Consider the set E of points a in B C (x 2 , R) such that a ≤ x 2 . Recall that a topological space B is connected i the only open and closed subsets A of B are A = ∅ and A = B. We use this denition with B = B C (x 2 , R) and A = E to exhibit a contradiction. The considered topology is now the induced topology on B C (x 2 , R). Making use of our assumption with the triple (a, b, c), where a

is in E, b = x 2 and c ∈ B(a, r) ∩ B C (x 2 , R), we see that for all a in E, B(a, r) ∩ B C (x 2 , R) ⊂ E. Since B(a, r) ∩ B C (x 2 , R) is a ball of B C (x 2 , R), E is an open subset of B C (x 2 , R). Again by our assumption, if a point c ∈ B C (x 2 , R) is not in E then B(c, r) ∩ B C (x 2 , R) cannot contain a point a ∈ E, otherwise the point c which is in 2.3. Existing total orders B(a, r), would be in E. It follows that B C (x 2 , R) \ E is also an open subset of B C (x 2 , R), hence E is an open and closed subset of B C (x 2 , R). Since B C (x 2 , R) is connected
, the open and closed subsets E and B C (x 2 , R) \ E cannot be both nonempty, a contradiction with the existence of the points x 1 and x 3 .

The order is going back and forth between point a, point b and point c. Note that in the particular case where X = R n>1 , this result is not a corollary of the fact that there is no continuous bijection between R n>1 and R. Indeed, not every order on R n>1 can be represented by a bijection on R (this would contradict the least-upper-bound property). This lemma tells us that for any total order in R n , functions supremum of two points, R n × R n → R n : (x, y) → sup(x, y) and inmum inf(x, y) present high irregularities with respect to the Euclidean metric, i.e. for any r, R > 0 there exists x 0 , x 1 , y 0 , y 1 ∈ R n such that ||(x 0 , y 0 ) -(x 1 , y 1 )|| < r and R < || sup(x 0 , y 0 ) -sup(x 1 , y 1 )||. Let I be an image composed of two objects of distinct colors with smooth boundaries. The rst object is composed of two values x 0 and x 1 with d(x 0 , x 1 ) < r. Points x 0 and x 1 are two versions of a color C x . The second object is composed by an unique value y with R < d(x 0 , y), representing color C y . Assume that sup(x 0 , y) = x 0 and sup(x 1 , y) = y. Regarding the object, the morphological dilation and erosion are inconsistent and might introduce aliasing on boundaries. The lemma is valid in any metric space and has strong negative implications. Namely, given a total order, it is always possible to nd an image where the erosion and dilation introduce important irregularities. An illustration of this phenomenon is given by the following toy example depicted in Fig. 2.2. The RGB image is composed of 3 dierent colors. Two close black (C x ) represented by a = (0, 0, 0) and b = (1, 0, 0), and a blue (C y ) represented by c = (0, 0, 255). According to the lexicographic ordering on coordinates, we are precisely in the situation described in the lemma. Fig. 2.2 shows us the result of a dilation using a 3 × 3 square as structuring element. Through the study of space lling curves, see [START_REF] Sagan | Space-lling curves[END_REF], the work [START_REF] Chanussot | Total Ordering Based on Space Filling Curves for Multivalued Morphology[END_REF] proposes several total orders on Z n that preserve as far as possible the notion of neighbourhood. For each point of the discretized multidimensional vector space, the spatial neighbourhood and the neighbourhood in the chain of the order are compared, as shown in Fig. 2.3 extracted from [START_REF] Chanussot | Total Ordering Based on Space Filling Curves for Multivalued Morphology[END_REF]. The average comparison between the spatial neighbourhood and the neighbourhood in the chain of the order gives a measure of the preservation of the neighbourhoods. This measure enables us to compare every total orders on Z n . [START_REF] Chanussot | Total Ordering Based on Space Filling Curves for Multivalued Morphology[END_REF] shows that the neighborhoods are signicantly better preserved by the Peano curve and the bit-mixing paradigm than by the lexicographic order. Despite the fact that the Peano curve gives slightly better results than the bit-mixing strategy paradigm, authors of [START_REF] Chanussot | Total Ordering Based on Space Filling Curves for Multivalued Morphology[END_REF] chose not to use the Peano curve because this order does not preserve the vectorial structure; i.e., in two dimensions, the point (0, 8) can be greater than the point (15,15).

In the present work we chose to focus the study only on topological considerations, without taking into account the algebraic structure of the value space. Note that the algebraic structure is often present by default without link with the physical quantity.

2.4. An order adapted to a given image

Motivation

As we have discussed above, it is not possible to create a total order that preserves neighbourhood on a multidimensional space. The philosophy of [START_REF] Chanussot | Total Ordering Based on Space Filling Curves for Multivalued Morphology[END_REF] and [START_REF] Flórez-Revuelta | Ordering of the RGB Space with a Growing Self-organizing Network. Application to Color Mathematical Morphology[END_REF] is to try to minimize the dierence between spatial neighbourhoods and neighbourhoods in the space of order. However it is possible to push this idea further. Even for the best total order in the sense of the measure dened in Fig. 2.3, our lemma tells us that the processing of a particular image can give highly irregular results.

As a consequence, it might be more interesting to look for the best order being given an image, instead of looking for the best order in general. Indeed, restricting the evaluation of a total order to a particular image, largely enhance the potential quality of the order.

An order on a multidimensional space can present important jumps that might not aect the processing of a particular image.

Firstly, given a specic image, the evaluation of a total order only gives importance to colors that lay in the image. Introduction of other (false) colors is due to the use of partial ordering. With total ordering, the values are chosen among the vectors of the input and as such there is no new colors. Indeed, no other colors are introduced by at erosion or dilation. Decreasing the size of the set to be ordered can lead to more regular orders. More precisely, it becomes easier to nd orders that avoid having two points close in the color space and far in the order. Secondly, problematic triplets of values described in the lemma can be tolerated to the extend that values do not occur in same zones of the image. Indeed, 2.4. An order adapted to a given image if two values never occur simultaneously in a structuring element, there relative position in the order will not aect the processing.

As explained in [START_REF] Chanussot | Total Ordering Based on Space Filling Curves for Multivalued Morphology[END_REF], it is clear that if points are close in the color space, they should remain close in the order. However, the reverse is not always required. Let us consider a binary image represented on the real line, where black and white are not represented by 0 and 1 but by 0 and 10. This situation presents two points close in the order chain and distant regarding the metric of the color space. However, this does not introduce any irregularity in the morphological operators. Thus the cost dened in Fig. 2.3 penalizes situations that are not problematic for mathematical morphology.

In order to transpose this topological intuition of closeness both in value space and in spatial space should imply closeness in the order, we introduce a cost function to be minimized by the total ordering.

Both works [START_REF] Ledda | Majority Ordering for Colour Mathematical Morphology[END_REF] and [START_REF] Lezoray | Graph-based ordering scheme for color image ltering[END_REF] proposed to restrict the order to values present in the image. Authors of [START_REF] Ledda | Majority Ordering for Colour Mathematical Morphology[END_REF] propose to order values according to their importance in the histogram of the image. While being a natural idea, it does not rely on any geometrical consideration on the value space. Authors of [START_REF] Lezoray | Graph-based ordering scheme for color image ltering[END_REF] propose to build local orders based on optimal Hamiltonian path.

They do not order the entire set of values of the image, and the associated dilation and erosion are not strictly speaking morphological operators. However the present work is in the continuation of their idea. We try here to build an order on the full value set present in the image, optimal with respect to a topological criteria.

Cost function

Given an image, we would like to dene a cost function that measures the quality of the order regarding rank based operators. We rst need to quantify the criteria of the previous section. We dene the following notion of co-occurrence of values a and b:

C I (a, b) = Card{p ∈ Ω | ∃q ∈ B(p), (I(p), I(q)) ∈ {(a, b), (b, a)}}.
(2.6)

The computation of co-occurrences involves that one has xed a typical size/shape of the structuring element B which will be used in subsequent processing. The motivation for introducing a notion of co-occurrence in the cost function is the following: if two colors are not present in the same zones of the image, they cannot produce aliasing. Let a and c be two close points in V according to a given metric distance d(a, c), such that a < c.

       a < c C I (a, b) and C I (c, b) signicant d(a, c) small d(a, b) signicant =⇒ a < c < b or b < a < c (2.7)
Adjectives signicant and small induce an increasing or decreasing behavior of the cost function with respect to the dierent quantities. Let us endow I(Ω) with a total order ≤. We can dene the following quantity:

P I (≤) = a, b, c ∈ I(Ω), a < c < b g(a, b, c) • (C I (c, a) ∧ C I (c, b)), (2.8)
where ∧ stands for min, with,

g(a, b, c) = f ((d(c, a) ∧ d(c, b)), d(a, b)) ,
and f (•, •) an increasing function according to the rst variable and decreasing according to the second.

Given an image I(p), this adapted cost function is more tolerant for some specic orders than the cost function dened in [START_REF] Chanussot | Total Ordering Based on Space Filling Curves for Multivalued Morphology[END_REF]. The cost function P I has been designed to represent as well as possible what is expected of an order. One of its main advantage is to take into account the above condition, see Eq. (2.7), weaker than the one required in [START_REF] Chanussot | Total Ordering Based on Space Filling Curves for Multivalued Morphology[END_REF], see Fig. 2.3. However, as a standard image often contains more than ten thousand dierent colors, this cost function presents the serious drawback of not being computable. Thus, given two orders, it is dicult to compare them using this cost function. Nevertheless, it is possible to try to minimize this cost function using a recursive procedure, without computing globally the cost for the full set of points.

Minimisation of the cost function

Overview of the algorithm

The idea is to divide the set I(Ω) into a collection C of clusters and to compute an optimal order on C considering each cluster as a single point. Then each cluster is ordered individually. The orders on individual clusters are merged according to the order on C to obtain an order on I(Ω). The point of the clustering is to make this operation possible by reducing the number of parameters of the minimization. Here are the main steps of the algorithm.

Perform a clustering of the data I(Ω) in a given number n of clusters: {Cluster i } 1≤i≤n .

Order the n clusters according to C(i, j), D(i, j), and the corresponding cost function P where D(i, j) represents the minimum distances between clusters i and j. C(i, j) represents the co-occurrence of clusters i and j in dierent structuring elements.

Perform the same procedure recursively on each cluster.

After a given threshold, stop the recursion. For each remaining cluster, order all its points according to a criterion based on distances to previous and next clusters.

Merge the orders: for x in Cluster i and y in Cluster j , x < y if and only if Cluster i < Cluster j , or i = j and x < y in Cluster i .

In this short description we did not mention a signicant source of complications. Indeed, when the recursive procedure is applied to Cluster i , one has to take into account neighbour clusters. Indeed, at the rst step of the recursion, nothing has to be taken into account except the considered set of colors I(Ω). However, unlike the set I(Ω), Cluster i can no longer be considered as isolated from the rest of the color values. If there exist colors c k in Cluster i and colors c l in Cluster j such that d(c k , c l ) is small, then it is not possible to order Cluster i without taking into account Cluster j . To order Cluster i , one needs to know the set of its neighbour clusters and their relative ordering. These are the main ingredients of the recursive algorithm. The algorithm uses ve functions.

A function Neighbour input: a list S 1 of clusters and a Cluster index in S 1 output: a sublist S neighbour of N neighb elements of S 1

The elements of S neighbour are selected from the list S 1 according to their co-occurrence with Cluster index . The co-occurrence between the clusters are computed from the image Ĩ where Ĩ(p) = i if and only if I(p) ∈ Cluster i . We impose Cluster index ∈ S neighbour . N neighb is an arbitrary parameter of the algorithm, see the result section.

A function IndexCutting input: a list S 1 of clusters and a Cluster index in S 1 output: a list S Cluster index of N Clust sub-clusters of Cluster index . N Clust is an other arbitrary constant of the algorithm, see the result section.

Call the function Neighbour to replace S 1 by the shorter list S neighbour . For each cluster Cluster i in S neighbour , nd the point c i in Cluster index that minimizes the distance to Cluster i . If i = index, let c i be the barycentre of Cluster index . The set {c i } is completed with random points in Cluster index to reach a minimum number of points. Perform a clustering of Cluster index using a k-means algorithm initialized with the {c i }. Return S Cluster index the new set of clusters.

A function MainOrder (To be used when the recursion depth is ≤ threshold) input: a list S 1 of clusters with an order < 1 , a Cluster index in S 1 and a list S Cluster index of sub-clusters of Cluster index output: an order < 2 , on the list S 2 = S 1 ∪ S Cluster index \ Cluster index .

Perform a minimization of P on S 2 , using the co-occurrences of clusters C(Cluster i , Cluster k ) and the distance D(Cluster i , Cluster k ), such that the new total order < 2 is compatible with the initial order < 1 . If Cluster i and Cluster j are in S 1 ,

(Cluster i < 1 Cluster j ) ⇒ (Cluster i < 2 Cluster j )
Furthermore, for Cluster i in S Cluster index and Cluster j in S 1 = S 1 \ Cluster index ,

Cluster i < 2 Cluster j ⇔ Cluster index < 1 Cluster j
The situation is summarized in Fig. 2.4. A function SimpleOrder (To be used when the recursion depth is > threshold) input: a list S 1 of clusters with an order and a Cluster index in S 1 output: an order on Cluster index . for all the sons in S index taken increasingly for < 2 do listorder ← Concat (listorder, FuncRecursive (S 2 , son, D + 1))

end end

The main program is simply the function FuncRecursive called with the list S 1 reduced to I(Ω) with Cluster index = I(Ω) and D = 0. The non-deterministic aspect of the clustering step induces that two successive minimizations of the cost function do not necessarily produce identical orders. Furthermore, the order dened on ε(I)(Ω) is not necessarily induced by the order dened on I(Ω). Thus, in order to preserve standards properties of mathematical morphology such as the idempotence of composed dilation and erosion, ε•δ, one has to use the same order for all operators on the complete processing chain.

Optimization over the permutation space

We note that, according to this algorithm the set to be ordered is no longer the set of all colors present in the image, but the set of clusters S at each level of recursion. The cost function P can be calculated if the cardinal n of S is reasonable. If the set S does not exceed n = 9 elements, it is conceivable to calculate the cost P of each permutation and select the minimum. However, this procedure requires an important computational time. An alternative solution consists in optimizing the cost function iteratively. It is well known that from a permutation it is possible to reach any other permutation by composing transpositions. Given a permutation, we select the transposition that minimizes the cost function, and we repeat the process until we fall in a minimum. As there is no proof that the cost function is convex, the minimum might only be a local minimum. However, under n = 9 elements it is possible to compare the result with the global exploration. In all tested situations with n = 9 elements, the optimization by transposition does not fall in a local minimum.

Computational cost

Given an order, the computation of the cost function P involves a sum over

n 3 = n(n -1)(n -2) 6 ,
terms where n is the number of points. If N Clust is the number of clusters produced by the function IndexCutting, the minimization of P is performed in (N Clust ) depth times. Our Matlab code launched on a standard computer gives 1, 25.10 -2 s for a single computation of P and 1.8 for its minimization over a set of 10 points (the second gure is an average over 50 runs). The choice of the clustering algorithm might induce signicant variation in the global computational cost. Using a k-means algorithm, N clust = 10, N neighb = 15 and depth = 2, our Matlab code requires 8 minutes to process an image of size 256 × 256 pixels composed of 40000 dierent colors.

Order invariances

The cost function P (I) does not depend on the coordinates of colors but only on their mutual distances and their co-occurrences. As the notion of co-occurrences remains unchanged under bijective transformations, the cost P (I) is invariant under any isometric transformation. Note that this is not the case for the lexicographic and the bit-mixing order. However, the choice of a particular function f (•, •) has an inuence on the class of invariances. For instance, the function

f (min i,j , d(c i , c j )) = min i,j d(c i , c j ) ,
(2.9)

where min i,j = d(c, c i ) ∧ d(c, c j ) also provides invariance of P (I) under homothetic transformations. We can also note that for the function

f (min i,j , d(c i , c j )) = min α i,j d(c i , c j ) β ,
(2.10) homothetic transformations simply result in the multiplication of P (I) by a positive constant. Consequently, as the notion of minimum is invariant under increasing transformation, the minimization of P (I) should remain relatively stable. If T is an isometric (or homothetic) transformation of the value space V, and φ a morphological operator {Ω, V} → {Ω, V}, then T and φ commute for any image I, i.e., φ(T (I)) = T (φ(I)).

Results of morphological image processing

For each image, the minimization is launched with the following parameters:

Function f : f (x, y) = x • G(y), (2.11)
where G(.) is a gate function with linear decrease, see Fig. 2.6. The parameter of the gate G is set to 1/10 of the diameter of the value space V.

The cost P becomes:

P = a < c < b (d(a, c) ∧ d(a, b)) × G(d(a, b)) × (C I (c, a) ∧ C I (c, b)).
(2.12)

At each level of the recursion, the set S is divided into n ≤ 10 clusters, assuming that its size allows it.

The number of neighbors selected by the function Neighbour is set to 15 at most, depending on the number of existing clusters. During the MainOrder function the cost P is evaluated on set of size 25 at most.

The recursion is stopped when depth = 2.

2D point sets

Before going into morphological operations on images, we present here a comparison of the lexicographic, the bit-mixing, and the adapted order on a few sets of points in a square of R 2 , without introducing the notion co-occurrences. Fig. 2.7 presents an example of ordering of 20 points, according to the lexicographic order, the bit-mixing order, and the adapted order. The adapted order is directly minimized over the set of permutations, without going through the clustering step.

Over a set of 200 points, the cost P can still be computed, which enables the comparison between the dierent orders. The inuence of the number of clusters, the depth of the recursion and the spatial conguration is studied in Fig. 2.8. Each experiment focuses on a specic geometric pattern. Points are randomly drawn around the geometric patterns.

The image adapted order is computed for a recursion depth varying from 1 to 4. At each depth, the adapted order and its cost are computed for a number of clusters varying on odd numbers from 2 to 12. For each conguration, the costs of the dierent total orders are averaged over 20 runs. The second column of Fig. 2.8 presents the costs depending on the depth. At each depth, the cost is averaged over the choice of the number of clusters. The third column represents the evolution of the cost according to the number of clusters at depth = 2. Fig. 2.8 enables to draw several conclusions. Except exceptional situations, the recursive minimization provides better results than the lexicographic and the bit-mixing orders. It succeeds in adapting to each situation. The decrease of the cost function with the depth lays mainly during the two rst steps of the recursion. This fact motivated the choice of depth = 2. Fig. 2.8 illustrates that the quality of the dierent orders strongly depends of the conguration of the point cloud. Indeed, although the lexicographic order seems to be the less interesting order, it is possible to nd congurations where it performs better than the bit-mixing and the adapted order. The behavior of the adapted order with respect to the number of clusters is more dicult to handle that the inuence of the depth. For rectilinear distributions, increasing the number of cuts of the point set seems to slightly decrease the quality of the adapted order. When the geometrical conguration has an intrinsic number of clusters, overestimating the number of clusters does not guarantee to obtain optimal results. In case of uniform distribution of the points, the cost seems to be globally decreasing with the number of clusters. The absence of assumption on the spatial conguration pushes to chose a signicant number of clusters. In the next experiments on images, 10 clusters has been judged as the best compromise between accuracy and computation time. In the dierent plots, the lexicographic order is in red, the bitmixing in green and the adapted order in blue.

Color Imaging

Unfortunately unlike the previous 2D point sets, the set of colors in a standard image is too large to enable the computation of the cost of a total order. It is then not possible to numerically evaluate the performance of the optimization. However, in the following examples the dependence on the order of the irregularities produced by the morphological processing is visually appreciable.

We present results of morphological color processing obtained for two dierent RGB images. The rst one is a microscopic blood vessel from a uorescence microscope, the second one is a natural color image. For both of them P is minimized by the recursive algorithm discussed in previous sections. The distance d(., .) between colors is the Euclidean distance of the RGB color space. give dierent results but are both coherent with the original image.

Using each total order, we can compute morphological color operators. Fig. 2.10 gives the corresponding openings and closings using as structuring element a square of 7 × 7 pixels. As expected, the lexicographic ordering produces important aliasing around the blue spots. The bit-mixing paradigm and the image adapted total order give dierent results but both preserve the regularity of boundaries. The same study is performed on the second example. Fig. 2.11 provides the original RGB image and the image representation of the three orders. Unlike the previous examples, lexicographic order is able to distinguish all the interesting objects of the image.

Furthermore, it seems to give an order smoother than the bit-miximg paradigm and the image adapted total order. As we can observe on Fig. 2.12, which gives the result of median ltering, the regularity of the grey-scale projection of the lexicographic order is only an illusion. Opening and closing operators of size 9 × 9 are depicted in Fig. 2.13 and Fig. 2.14. Both lexicographic and bit-mixing order present aliasing on blue and yellow boundaries. The adapted order presents a slight aliasing in the background, which is not present with lexicographic and bit-mixing orders. However on this example, our image adapted total order is the only one of the three orders that provides satisfying results in terms of regularity for opening and closing.

Diusion Tensor Imaging (DTI)

We recall that in DTI, each pixel of the image contains a 3 × 3 symmetric positive denite matrix, that is a point in the space SPD(3). A matrix of SPD(3) can be represented as an ellipsoid, thus DTI image can be visualized using ellipsoids, see the example in Fig. 2.15. However, this new framework is also based on partial order relationships. Similarly to the color section, we restrict the comparison of the irregularities to total orders.

As the matrix is symmetric, it can be parametrized by the 6 upper-diagonal coecients, see Fig. This parametrization of SPD(3) induces a lexicographic and a bit-mixing order, named respectively LEX1 and BM IX1 in what follows. Another approach consists in representing SPD(n) matrices with rotations and eigenvalues. A symmetric matrix can be diagonalized in an orthonormal basis. Then each SPD(3) matrix can be represented by 3 eigenvalues and a rotation matrix. Using any angular representation of the rotation matrix, the SPD(3) matrix can be represented by 3 eigenvalues and 3 angles. Let us consider such a parametrization where the eigenvalues are sorted decreasingly. The associated lexicographic order is named LEX2 and the bit-mixing order BM IX2.

Using our framework, the image adapted total order is calculated according to two metrics: the metric associated with the Frobenius scalar product, i.e., A, B = tr(AB t ); In every processing example, the structuring element B is a square of 5 × 5 pixels. Using the openings from Fig. 2.17 as markers, the geodesic reconstruction [START_REF] Soille | Morphological Image Analysis[END_REF] has also been computed. The results are depicted in Fig. 2.19 and as expected, this operator based on geodesic dilations will recover the contours of the objects which have not been suppressed by the opening.

Finally, the comparative result of median ltering is presented in Fig. 2.20.

Visually, the image adapted total orders give results at least as good as other orders, and better ones in several situations. The choice of the Euclidean or log-Euclidean metric does not seem to introduce important changes.

Total order in the Poincaré upper half plane

As discussed in previous sections, the proposed order is based only on the notion of distance, independently from the algebraic structure. This framework suits perfectly to the case of images whose pixels lies in a Riemannian manifold, where the distance is known. We present here a situation where the image is valued in such a manifold.

Mathematical morphology for images valued in uni-dimensional Gaussian laws has already been studied in [START_REF] Angulo | Morphological Processing of Univariate Gaussian Distribution-Valued Images Based on Poincaré Upper-Half Plane Representation[END_REF]. One of the most common distances on Gaussians is the distance induced by the Riemannian metric called the Fisher metric. This is a simple example where the value space is not a vector space. In the case of one dimensional Gaussian parametrized by their means µ and their standard deviation σ, the distance between two Gaussians G 1 and G 2 is given by:

d(G 1 , G 2 ) = cosh -1 (1 + 1 2 (µ 1 -µ 2 ) 2 + (σ 1 -σ 2 ) 2 2σ 1 σ 2 ).
In the mean/standard deviation half-plane, the shortest paths of the Fisher metric are half ellipses centered on the µ-axis. The mean/standard deviation half-plane endowed with the Fisher metric is a Poincaré upper half plane of curvature -1 2 .

The studied example is a time lapse sequence of grey-scale images from a retina. At each pixel, we dispose of 20 successive acquisitions. By assuming a Gaussian distribution on the successive acquisition, we obtain a Gaussian valued image represented in Fig. 2.21.

The observation of the distribution of Gaussians in the upper half plane presented in Fig. 2.21(c) strongly suggests to endow the space of Gaussian laws with the Fisher metric.

The proposed solution to obtain total orders is currently the only total order framework that takes the geometry of the space into account.

Conclusions and Perspectives

To our knowledge, this work is the rst study that rigorously formulates and demonstrates the irregularity issue of total orders in metric spaces. Given a total order, the introduced Lemma tells us there are always images where morphological operators introduce irregularities and aliasing. Given an image and a total order, we exhibit the problematic triplets of values. The identication of problematic triplets leads us to a cost function measuring the quality of a total order regarding the back and forth problem. Due to the restriction of the problem to a specic image and to a precise identication of the problematic situations, we know that the order minimizing the cost function produces less irregularities than the existing propositions, namely [START_REF] Chanussot | Total Ordering Based on Space Filling Curves for Multivalued Morphology[END_REF][START_REF] Flórez-Revuelta | Ordering of the RGB Space with a Growing Self-organizing Network. Application to Color Mathematical Morphology[END_REF].

As the cost function is not explicitly computable due to the size of the set to be ordered, we propose a recursive minimization procedure based on successive clustering of the set to order, so as to nd a total order adapted to the image. Future work might be led on boundaries regularity measurement so as to provide a quantitative measure of the quality of the minimization. A reasonable lead would be to evaluate the irregularity of the level-sets channel-wise and merge the results.

A strength of the proposed framework is its generality. Indeed, for a large majority of images, the value space is equipped with a metric. The image adapted order can thus be computed for almost any image. Note that for a given type of images, the adapted metric is often open to debate, see for instance [START_REF] Ledoux | How to specify or identify the most accurate multispectral distance function for mathematical morphology?[END_REF].

We have shown the interest of our method on some examples, however the bit-mixing total order proposed in [START_REF] Chanussot | Total Ordering Based on Space Filling Curves for Multivalued Morphology[END_REF] remains an interesting solution in the Euclidean case. Indeed our current minimization procedure does not yet provide average results signicantly better than the bit-mixing, while the latter is independent from the image and requires no pre-processing.

The minimization of the cost function is subject to several potential improvements, on its accuracy on the one hand, and on the computation speed on the other hand. The use of techniques such as evolutionary algorithms might for instance enable to signicantly enhance the number of clusters at each step of the recursion, see [START_REF] Flórez-Revuelta | Topology-Preserving Ordering of the RGB Space with an Evolutionary Algorithm[END_REF].

Thus our future research will be essentially focused on the minimization of the cost function.

N-ary morphology

Abstract

Mathematical morphology on binary images can be fully described by set theory. However, it is not sucient to formulate mathematical morphology for grey scale images. This type of images requires the introduction of the notion of partial order on grey levels, together with the denition of sup and inf operators. More generally, mathematical morphology is now described within the context of the lattice theory. For a few decades, attempts are made to use mathematical morphology on multivariate images, such as color images, mainly based on the notion of vector order. However, none of these attempts has given fully satisfying results. Instead of aiming directly at the multivariate case we propose rst an extension of binary mathematical morphology to an intermediary situation: images composed of a nite number of independent unordered labels. We propose then a second extension to a continuous case.

This chapter is mainly adapted from [START_REF] Chevallier | Probability density estimation on the hyperbolic space applied to radar processing[END_REF].

Introduction

The philosophy of the previous chapter was to adapt the data to the theory of mathematical morphology. The aim was to put an interesting structure on data, a total order in this case, that would enable the use of morphological operators. However in the multivariate case, addressing data only through an order, partial or total, is often not fully compatible with the original structure. Product orders for instance are equivalent to marginal processing. Lemma 2.3.1 shows that total orders on R n>1 are never compatible with distance information in morphological processing.

Instead of adapting the structure of the data to the mathematical morphology theory, one could try to adapt the theory to the data. This chapter do not aim directly at the full multivariate case. We propose to adapt the morphological framework rstly to images valued in n independent categories and secondly to images valued in mixtures of this categories.

Let us start by quoting [START_REF] Serra | Image analysis and mathematical morphology[END_REF]. Objects in space generally have three dimensions, which are reduced to two dimensions in a photograph or on the retina. In this projection the luminance of the point located along a line oriented directly away from the viewer are not summed, because most physical objects are not translucent to light rays, in the way they would be to X ray but are opaque. Consequently, any object that is seen hides those that are placed beyond it with respect to the viewer: this self-evident property is a basic one. In fact it serves as a starting point for mathematical morphology, since, whenever we wish to describe quantitatively phenomena in this domain, a set-theoretic approach must be used.. Here is one of the main dierences between (nonlinear) image processing and sound processing. In sound processing, audio sources of the same environment are summed, whereas opaque objects of a 3D scene hide themselves. In image processing we might not want to mix colors of neighbor objects, as does the linear convolution, because this mix might have no meaning in terms of objects. Thus this key idea can be reformulated by saying that mathematical morphology is based on the reduction and the extension of the surface of the dierent objects over their neighbors. That leads naturally to the two basic morphological operators in binary images.

Binary images. In such images, there are only two kinds of objects: black or white objects. Two dual and adjoint operators have been dened: the erosion and the dilation.

The erosion extends the black objects over the white objects, the dilation extend the white objects over the black objects. Formally, a binary image can be seen as a support set Ω, and X a subset of Ω. Set X represents the white points of the image and the complement X = X c the black ones. We assume here that Ω is a subset of R n , which is generally the case. Let B be a subset of Ω called the structuring element. We assume that Ω has a translation operation. The erosion ε B (X) and the dilation δ B (X) of X according to a structuring element are dened as follows, see [START_REF] Matheron | Random sets and integral geometry[END_REF]:

ε B (X) = y∈B X -y = {p ∈ Ω : B p ⊂ X} = x : ∀p ∈ B, x ∈ X p , (3.1) δ B (X) = y∈B X y = {x + y : x ∈ X, y ∈ B} = p ∈ Ω : X ∩ Bp = ∅ , (3.2)
where X = {-x : x ∈ X} is the transpose of X (or symmetrical set with respect to the origin O) and X p = {x + p : x ∈ X} the translate of X by p. We note that erosion and dilation correspond respectively to the Minkowski substraction and addition of set X by B; i.e., ε B (X) = X B and δ B (X) = X ⊕ B. For the sake of simplicity, we limit the rest of our notation to symmetric structuring elements: B = B.

Grey-scale images. With the appearance of grey-scale images, mathematical morphology was reformulated in terms of inf and sup convolution where the kernel is the structuring element B, see [START_REF] Serra | Image analysis and mathematical morphology[END_REF]. That corresponds to the so-called at case, where the kernel is just a shape. It is also possible to use structuring functions, but this is out of the scope of this work. An image is now considered as a function I dened as

I : Ω → V p → I(p)
where V is the set of grey-levels, which can be generally assumed as a subset of the real line V ⊂ R. Grey-scale at erosion and dilation of I by structuring element B are now dened as follows:

ε B (I)(p) = inf q∈Bp {f (q)} , (3.3) δ B (I)(p) = sup q∈Bp {f (q)} . (3.4)
In this classical framework, each grey-level is not fully considered as an independent label (i.e., a dierent category) but simply as an intermediary level between black and white.

This point of view is actually justied when interesting objects of the images are local extrema. A mathematical way to see this formalism is to note that at erosion and dilation commute with processing upper level set by upper level set, that is operators (3.3) and (3.6)

are respectively equivalent to while the erosion restores the dark structures which are still present.

ε B (I)(p) = sup h ∈ R : p ∈ X + h (I) B , (3.5) δ B (I)(p) = sup h ∈ R : p ∈ X + h (I) ⊕ B ,
In the current example, it is not possible to remove the central grey spot using erosion and dilation with B larger than the spot size. This grey spot is not considered as an interesting object in itself but simply as an intermediary value between the black object and the white object. If this assumption is often coherent, this is not always the case.

Let consider the grey-scale image in Fig. 3.2(a). In this image, each grey level has the same semantic level: each represents a dierent component, sometimes called a phase.

However, in the morphological grey-scale framework, the grey is processed as an intermediary level. It is possible to replace each grey level by a color, see Fig. 3.2(b). We would like then to process both images using the same approach. Multivariate images. For multivariate images, no canonical framework has yet appeared. Most processing consist in endowing the structure with a partial order relationship.

The structure has to be a complete lattice in order to dene erosion and dilation in terms of inf and sup. The notion of order induces the notion of intermediary level as in the precedent framework. It is rather intuitive to consider grey as an intermediary between black and white. The situation is dierent for colors: using the product ordering in the RGB space red is an intermediary color between black and white. This assumption often has no meaning in term of image processing. In an image, the color red usually has the same semantic level as black or white: the color red represents a class of objects, it is not a transition between black objects and white objects. If the phenomenon already exists in the grey-scale framework, it turns out now to be a signicant problem, see example in Aim and chapter organization. The historical way of extending the mathematical morphology theory from binary to multivariate images is summarized in Fig. 3.3(a).

However, the gap between grey-scale and multivariate is much more signicant than the gap between binary and grey-scale. As discussed previously, the grey-scale structure only enables to dene intermediary levels between two references. This structure is obviously too weak to describe a multivariate information. Avoiding the notion of order becomes one of the goals of recent developments of mathematical morphology. Authors of [START_REF] Gronde | Sponges for Generalized Morphology[END_REF] propose to lighten the lattice structure, while authors of Carlinet and Géraud (2015) build a tree of shapes on multivariate images without requiring a total orders. Although very dierent from these propositions, our work is based on the same underlying idea.

Before extending mathematical morphology to multivariate images, we might rst want to dene a coherent approach for mathematical morphology with n independent unordered labels, without considering them as intermediary levels. This is the aim of the rst part of the chapter. The dierence between the frameworks can be interpreted in term of a change of metric on the value space:

grey-scale framework:

∀(i, j), d(label i , label j ) = |i -j|; n-ary framework: ∀(i, j), d(label i , label j ) = 1.
This work is not the rst to consider the problems of classical mathematical morphology for images composed of independents labels. Authors of [START_REF] Busch | Morphological Operations for Color-Coded Images[END_REF] and [START_REF] Ronse | Morphology on Label Images: Flat-Type Operators and Connections[END_REF] have very similar motivations but the development we propose is dierent. In contrast to operators proposed in [START_REF] Busch | Morphological Operations for Color-Coded Images[END_REF], [START_REF] Ronse | Morphology on Label Images: Flat-Type Operators and Connections[END_REF] or labeled openings from [START_REF] Hanbury | Morphological operators on the unit circle[END_REF], we are interested in lling gaps left by anti-extensive operators. We note that the theory of morphological operators for partitions introduced in Ronse ( 2013) and hierarchies of partitions in [START_REF] Meyer | Adjunctions on the lattice of hierarchies[END_REF] are not compatible with our framework.

The second part of the chapter aims to extend the morphology dened in the rst part on independent (discrete) labels to a continuous case. This part of the theory is less elaborated than the discrete case but provides already the key ingredients for future developments. The continuous setting involves to work with images valued in the (n -1)simplex. The approach adopted here will be based on a fuzzy representation and limited initially to the edges of the simplex, i.e., only two labels (or categories) contributes at each pixel. We should point out that our fuzzy framework is not related to the theory of fuzzy morphology as developed for instance by [START_REF] Deng | Grey-Scale Morphology Based on Fuzzy Logic[END_REF]; [START_REF] Bloch | Duality vs. adjunction for fuzzy mathematical morphology and general form of fuzzy erosions and dilations[END_REF].

In parallel to the approach considered here for the (n -1)-simplex, another recent work by [START_REF] Franchi | Ordering on the Probability Simplex of Endmembers for Hyperspectral Morphological Image Processing[END_REF] has considered the extension of morphological operators using partial orders on this space.

The rest of the chapter is organized as follows.

A proposition of n-ary morphological operators and a study of their theoretical properties is done in Sections 3.2 and 3.3. Some applications to image ltering are discussed in Section 3.4.

Section 3.5 proposes a continuous extension of n-ary morphology for images arising from fuzzy classications.

Section 3.6 concludes and closes the chapter.

n-ary morphological operators

Let us come back to the key idea of mathematical morphology is to reduce and extend objects over their neighbors. In the case of binary images, two operations were introduced:

the erosion extends the black over the white and the dilation extends the white over the black. In a general way, we would like to allow to reduce and extend the surface of each label of object. This makes four theoretic operations in the binary case, reduced to two in practice due to the coincidence of certain operations: reducing the black is the same as extending the white and conversely. This duality is one of the basic principle of binary morphology.

Dilation and erosion of label i

Let I be an n-ary image dened as

I : Ω → {1, 2, • • • , n} p → I(p)
In the n-ary case, it seems natural to try to introduce the corresponding pair (ε i , δ i ) of operators for each label i. Erosion ε i is the operator that reduces the surface of the objets of label i, and dilation δ i the operator that extends the label i. Above n > 2, we unfortunately lose the duality between operations, such that the number of elementary operators is then equal to 2n. Let us formulate more precisely these operators.

The dilation of label i on image I by structuring element B presents no diculty: Let us rst dene the erosion for pixels where there are no ambiguities. Thus the erosion of label i on image I by structuring element B is given by

δ i (I; B)(x) = I(x) if ∀p ∈ B x , I(p) = i i if ∃p ∈ B x , I(p) = i
ε i (I; B)(x) =    I(x) if I(x) = i i if ∀p ∈ B x , I(p) = i θ(x, I) otherwise (3.8)
We will address later denition of θ(x, I). Sections 2.2 and 2.3 are independent of θ. Although the image is a partition of Ω the proposed framework diers from [START_REF] Ronse | Morphology on Label Images: Flat-Type Operators and Connections[END_REF] and Ronse (2013).

Opening and closing of label i

Once the dilation and erosion have been dened, we can introduce by composition of these two operators the opening and the closing on I by B of label i respectively as

γ i (I; B) = δ i • ε i = δ i (ε i (I; B); B) , (3.9) ϕ i (I; B) = ε i • δ i = ε i (δ i (I; B); B) .
E I i = I -1 (i). To simplify, 1 E I i will be noted 1 I i .
Due to the fact that operators (3.7) and (3.8) are not adjoint, classical theory of opening and closing by adjunction, see [START_REF] Serra | Image Analysis and Mathematical Morphology[END_REF]; [START_REF] Heijmans | Morphological image operators[END_REF], does not apply to our framework. Hence, we need to explore which properties still hold and in particular, if operators (3.9) and (3.10) are idempotent. Indeed, we have such a property of stability.

Proposition 3.2.1 Opening and closing of label i are idempotent operators, i.e.,

γ i (I; B) = γ 2 i (I; B), ϕ i (I; B) = ϕ 2 i (I; B).
Proof Since the binary opening is idempotent, one has

E γ i (I) i = E γ 2 i (I) i . Furthermore we have that E γ i (I) j ⊂ E γ 2 i (I) j
, for all j = i. Since sets (E i ) i from a partition of the support space, necessarily E

γ i (I) j = E γ 2 i (I) j
, ∀j. Indeed, if all the elements of a partition are extensive, then they all remain stable. Then γ i = γ 2 i .

Properties cannot directly be transported by duality, as in binary morphology, however the property remains true for the closing. We rst show the binary property εδε = ε. The binary erosion and opening can be written as

ε B (X) = ∪ Bx⊂X {x}, and γ B (X) = ∪ Bx⊂X B x . Then ε(γ(X)) = ∪ Bx⊂γ(X) {x}. Since {B x ⊂ γ(X)} = {B x ⊂ ∪ Bx⊂X B x } = {B x ⊂ X}, then ε(γ(X)) = ε(X). Thus, εδε = ε and by duality, δεδ = δ. Then E δ i i = E δ i ε i δ i i . It can be shown that E δ i j ⊂ E δ i ε i δ i j
for all j = i. Using the same reasoning as in the proof for the opening, we have that for all j, E δ

i j = E δ i ε i δ i j . In other words, δ i ε i δ i = δ i . Thus ε i δ i ε i δ i = ε i δ i , or equivalently ϕ i = ϕ 2 i .

Composed n-ary lters

We can now try to dene label lters from the openings and the closings of label i. In binary morphology, the simplest lters are of the following form: γ • φ and φ • γ. In the n-ary framework, with n = 2, they can be rewritten as The opening removes peaks smaller than the structuring element and the closing removes holes, which are dual notions in binary morphology. However, peaks and holes are no longer a dual notion in n-ary morphology with n > 2. Fig. 3.5 illustrates the dierence between openings and closings on a 3-ary image: three colors, black, white and red. The structuring element is a square whose size is half of the width of the red line. Removing the red line using φ black requires a structuring element twice bigger than with γ red . As a good candidate to lter out small object of a labeled image I, independently of the label of the objects, we introduce the operator ψ, named composed n-ary lter by structuring element B, dened as

γ 1 • γ 2 = φ 2 • φ 1 , and γ 2 • γ 1 = φ 1 • φ 2 .
ψ(I; B) = γ n (I; B) • γ n-1 (I; B) • • • • • γ 1 (I; B). (3.11)
Unfortunately on the contrary to γ •φ in binary morphology, ψ is generally not idempotent. Worst, the sequence ψ k do not necessarily converge. However we still have a stability property for relevant objects. Let us be more precise.

Proposition 3.2.2 Let Ω be a nite set. Given a structuring element B, the interior with respect to B of the composed n-ary lter ψ(I; B) converges for any image I, i.e., ∀i, ε(E

ψ k i ) converges . Proof Since ε = ε • δε, ∀i, ε(E ψ k i ) = ε(E γ i •ψ k i ). Furthermore, since ε(E ψ k i ) ⊂ ε(E γ j •ψ k i ),
we have that ∀i, ε(E

ψ k i ) ⊂ ε(E ψ k+1 i ). Since Ω is a nite set, ε(E ψ k i ) converges.
This property ensures that the variations between ψ k and ψ k+1 do not aect the interior of objects and is only limited to boundaries. Nevertheless, as we shown in section 4, ψ k is almost always stable after a few iterations.

n-ary geodesic reconstruction

The binary reconstruction can be transposed in the n-ary framework as follows. The proposition of reconstruction is similar to the one proposed in [START_REF] Ronse | Morphology on Label Images: Flat-Type Operators and Connections[END_REF].

Given two labeled images R and M , for each label i, Perform a binarisation of the reference R and the marker M between i and i, which correspond respectively to binary images X i and Y i . Compute γ rec (X i ; Y i ), that is the binary geodesic reconstruction of the marker in the reference.

Then, the n-ary geodesic reconstruction of the reference R by the marker M is given by The aim of this denition is to symmetrize labels. In Fig. 3.6 (d), the grey object is considered as an object in itself. The proposed reconstruction is a connected operator in the sense of [START_REF] Salembier | Flat zones ltering, connected operators, and lters by reconstruction[END_REF] 3.2.5. n-ary morphological gradient

γ rec (R; M )(x) = i if x ∈ γ rec (X i ; Y i ) M (x) if ∀i, x / ∈ γ rec (X i ; Y i )
The classical symmetric morphological gradient of a binary or grey-scale image I is dened as the dierence between the dilation and the erosion:

grad(I)(x) = δ B (I)(x) -ε B (I)(x),
where B is typically the unit ball. This concept can be extended to an n-ary image I by the following formulation of the n-ary morphological gradient:

grad(I)(x) = 1 δ 1 (I;B) (x), 1 δ 2 (I;B) (x), • • • , 1 δn(I;B) (x) .
(3.12)

where

1 δ i (I;B) (x) stands for 1 E δ i (I;B) i (x) and grad(I)(x) is an element of {0, 1} n . As
the standard gradient, a large part of its information is contained in its norm. Here,

grad(I)(x) 1 = i 1 δ i (I) (x)
tells us how many dierent labels lays in the neighbourhood of each point x.

Fig. 3.7 shows the dierences between the classical morphological gradient and the n-ary morphological gradient of a labeled image. We note that for the classical color gradient, the involved dilation/erosion are computed using a lexicographic (total) ordering. Both gradients have similar aspect. However, their respective grey levels have dierent meanings.

In the classical framework, the intensity represents a dierence of position in the order, whereas in the n-ary framework, the intensity represents the local variety of labels.

On the choice of an erosion of label i

Before any application, we need to come back to the erosion problem. More precisely, we need to dene a consistent rule to ll the space created by the erosion operation.

First of all, we note that the denition in Eq. (3.8) of the erosion ε i of label i does not indicate how to behave on the following set:

A = {x | I(x) = i and ∃p ∈ B x such that I(p) = i} .
For points x ∈ A we have to decide by which label to replace label i and therefore to dene ε i on A, i.e.,

ε i (I; B)(x) =    I(x) if I(x) = i i if ∀p ∈ B x , I(p) = i ? if x ∈ A (3.13)
Many alternatives are possible. Two criteria have to be taken into account: (i) the direct coherence in terms of image processing, and (ii) the number of morphological properties veried by the erosion, such as ε i (I; kB) = ε k i (I; B) where kB = {kx | x ∈ B} (i.e., homothetic of size k). Let us consider in particular the three following rules for x ∈ A:

1. Fixed-label erosion: Erosion always lls the gaps with label 1 (or any other xed label):

ε i (I; B)(x) = 1.
ε i (I; B)(x) = min(arg min j =i d j x ).
(3.16) Let us formalize the iterative behavior of the distance-based erosion by the following result on isotropic structuring elements.

where d j x = inf { x -p Ω | p ∈ Ω, I(p) = j}.
Proposition 3.3.1 Let (Ω, d) be a compact geodesic space and I a n-ary image on Ω. For any R 1 , R 2 > 0 the distance-based erosion of label i satises:

ε i (I; B R 1 +R 2 ) = ε i (I; B R 2 ) • ε i (I; B R 1 ).
where B R is the open ball of radius R.

Proof For the sake of notation, let X = E I i and X j = E I j , for j = i.

Let X = {x ∈ X|d(x, X) ≥ R 1 }, X is the binary eroded of X. Let proj X (a) = {b|b ∈ X, d(a, b) = d(a, X)} with A the closure of A. Let P = {X j } and for a ∈ Ω, I P a = {j|∃b ∈ proj X (a) ∩ X j }. Let X j = {a / ∈ X |proj X (a) ∩ X j = ∅} and P = {X j }.
Using the property of the min function, min(A ∪ B) = min({min(A), min(B)}), it can be shown that I P a = I P a for all a ∈ X implies proposition 3.3.1.

Let j ∈ I P a . Thus, ∃b ∈ proj X (a) ∩ X j . Let γ be a geodesic joining a and b with γ(0

) = b and γ(d(a, b)) = a. Let c = γ(R 1 ). For all x ∈ γ([0, R 1 [), x /
∈ X and b ∈ proj X (x). Hence x ∈ X j and c ∈ X j . The triangle inequality gives c ∈ proj X (a).

Indeed, assuming that c / ∈ proj X (a) easily leads to b / ∈ b ∈ proj X (x). Thus, j ∈ I P a .

Let j ∈ I P a . Thus, ∃c ∈ proj X (a) ∩ X j . By denition of the closure,

∃c n → c, c n ∈ X j . c n ∈ X j ⇒ (∃d n ∈ proj X (c n ), d n ∈ X j ). The compact assumption enables us to consider that d n → d (it is at least valid for a sub-sequence). d ∈ X j . By continuity, d(d, c) = lim(d n , c n ) = lim( X, c n ) = d( X, c). Thus, d ∈ proj X (c). Proposition 3.3.2 tells us that d ∈ proj X (a). Hence j ∈ I P a .
Thus, I P a = I P a .

3.4. Applications to image ltering Proposition 3.3.2 Using the notation introduced in the demonstration of proposition 3.3.1:

∀a ∈ X , ∀c ∈ proj X (a), ∀d ∈ proj X (c), d ∈ proj X (a). Proof Let a ∈ X , c ∈ proj X (a), d ∈ proj X (c). Since X is closed, X is open and the existence of geodesics implies that c / ∈ X . Hence c ∈ X and d(c, X) ≥ R 1 . Since c is a projection on X , we have c ∈ X , d(c, X) ≤ R 1 . Hence d(c, X) = R 1 and d(c, d) = R 1 . Let b ∈ proj X (a) and γ a geodesic such that γ(0) = a and γ(d(a, b)) = b. Let c = γ(sup{t|γ(t) ∈ X }). Since X is closed, c ∈ X . We have d(c , X) ≥ R 1 . Since c ∈ X , d(a, c ) ≥ d(a, X ) = d(a, c). Thus, d(a, b) = d(a, c ) + d(c , b) ≥ d(a, c) + d(c , X) = d(a, c) + d(c , X) ≥ d(a, c) + R 1 = d(a, c) + d(c, d) ≥ d(a, d). Hence d ∈ proj X (a).
Note that in a compact subset of a vector space, since proposition 3.3.1 is valid for any norm, the property holds for any convex structuring element. Note also that result 3.3.1 is based on a metric denition of the erosion. For vector spaces norms, this erosion is identical to the translation based erosion. Proposition 3.3.1 strongly links the distancebased erosion to a wave-front propagation. Without entering the details, we present the general idea of wave-front propagation and the relation with n-ary morphology. Let the curve C 0 (s) = (x 0 (s), y 0 (s)) be the contour of a planar shape S. Let C(s, λ) be the curve dened by the following partial dierential equation:

∂C(s,t) ∂t = α - → N (s, t) C(s, 0) = C 0 (s) (3.17)
where N is the normal vector at C(s, t). Authors of [START_REF] Arehart | Mathematical morphology: The Hamilton-Jacobi connection[END_REF] and [START_REF] Sapiro | Implementing continuous-scale morphology via curve evolution[END_REF] point out that C(s, λ) corresponds to the dilated or the eroded, depending on the choice α = 1 or α = -1, of S by a ball of radius λ. The underlying idea is that a dilation (resp. erosion) can be decomposed into a succession of innitesimal dilations (resp. erosion). The situation is slightly more complicated in the n-ary case, see Fig. 3.9. The propagation of each labels j = i in the gap left by label i during ε i is up to a certain point described by a proper wave-front propagation following Eq.(3.17), see Fig. 3.9(b). However note that in Fig. 3.9(b), if Front 1 follows Eq.(3.17), it is not the case of Front 2. Front 2 is elongated following the trace left by the point p.

According to this discussion, in all what follows we adopt the distance-based erosion.

Applications to image ltering

In the rst case study, depicted in Fig. 3.10 and Fig. 3.11, we consider the behavior and interest of the composed n-ary lter ψ(I; B). The aim is to lter out objects smaller than the structuring element B of the 4-ary labeled image of label i. On the one hand, note that in all the cases, the iterated lter converges rather fast to a stable (idempotent) result and that the dierence between the rst iteration and the nal result are rather similar. On the other hand, the dierent permutations produce dierent results, however in all the cases, the small objects seems better removed than in the case of the classical color order operators.

The second example, given in Fig. 3.12 attempts to regularize the 3-ary image, by removing small objects without deforming the contours of the remaining objects. More precisely, Fig. 3.12(a) represents the electron microscopy image of a ceramic eutectic, with three dierent phases after segmentation. The ltering process is composed of two steps: morphological size lter followed by geodesic reconstruction. We compare the results of ltering of the color image according to two pipelines: (i) color total order framework, Fig. 3.12(b), where the lter is an opening by reconstruction composed with a closing by reconstruction; (ii) n-ary framework, Fig. 3.12(c), where the marker is a n-ary lter ψ(I; B) followed by a n-ary geodesic reconstruction. In the case of the color ordering, black and blue are extreme color whereas red is the intermediary color. As we can observe, using the order-based approach all red objects that lay between black and blue objects are not extracted. Both n-ary and color total order frameworks give the same results for the blue grains. This corresponds to the fact that, in the color order, blue is an extreme color whereas red is an intermediary. Therefore the n-ary framework provides a more symmetric processing of all the labels.

The last example is a classication image from the brain.

; B) = γ 1 • γ 2 • γ 3 • γ 4 , (b) ψ(I; B) = γ 2 • γ 1 • γ 3 • γ 4 , (c) ψ(I; B) = γ 3 • γ 4 • γ 1 • γ 2 , (d) ψ(I; B) = γ 4 • γ 3 • γ 2 • γ 1 .
a classication where the red represents the grey matter, the green represents the white matter and the blue represents the cerebrospinal uid. The processing is same as for the second example. The miss-classied white matter around the brain and some missclassied grey matter spots around the cerebrospinal uid are successfully removed by the n-ary framework, whereas they remain after the classical processing.

A rst extension to a fuzzy case

A natural candidate space to generalize n-ary morphology is the (n -1)-simplex, that is to say the space of discrete probability distributions on n elements. Instead of assigning one and only one label per pixel, we assign now a discrete probability distribution representing a local mixture of classes.

In this setting, an image is a function of the form:

I : Ω → {(a 1 , • • • , a n ), a i ≥ 0, i a i = 1} p → I(p)
The aim of this section is to adapt to the continuous n-ary setting, the denitions of the dilation and erosion per label as in (3.7) and (3.8). As we have pointed out in the introduction, the approach developed here has no direct link with the classical grey-scale fuzzy morphology, see for instance [START_REF] Deng | Grey-Scale Morphology Based on Fuzzy Logic[END_REF]; [START_REF] Bloch | Duality vs. adjunction for fuzzy mathematical morphology and general form of fuzzy erosions and dilations[END_REF] as relevant works on the latter topic.

Basic operators and their properties

Similarly to (discrete) n-ary case the dilation is more straightforward than the erosion. Let I(p) i be the i-th element of I(p), and I i be the real valued image associated to the i-th component. A natural way to dilate the component i is to perform a grey-scale dilation of the component i, followed by a rescaling of the other components in order to belong to the (n -1)-simplex: with

δ i (I; B)(x) j =      δ(I i ; B)(x) if j = i 1-δ i (I;B)(x) i 1-I(x) i .I(x) j if j = i (3.18) (a) (b) (c) (d) (e) (f ) (g) (h) (i)
1 -δ i (I; B)(x) i 1 -I(x) i = 1, when δ i (I; B)(x) i = I(x) i = 1.
One could be tempted to dene the erosion the same way, that is to erode the component i and rescaling the other components. In the n-ary case this is equivalent to leaving the gap of Fig. 3.4(b) empty. The surface of the label i is reduced, but not replaced by a specic label, which is not entirely satisfying neither consistent with the discrete case.

Furthermore, contrary to the dilation the rescaling is impossible when I(x) j = 0 for all j = i.

The amount to be lled is I(x) i -ε(I i ; B)(x), where ε(I i ; B) stands for the classical erosion of the partial domain dened function I i . One has to decide with which label or set of labels to complete the distribution. Unlike in section 3.2.2, the stability properties of composed operators of the form ϕ i = ε i • δ i or γ i = δ i • ε i strongly depends on that choice.

Following the philosophy from Section 3.3, the components used to ll the gaps should be chosen among components present in the local structuring element B x . The diculty is that many of the possible choices do not provide idempotence of ϕ i and γ i .

Consider for instance the following erosion rule. A grey-scale erosion is performed on component i. At x, let j = i be the index of the component with the highest value in B x . When j is not uniquely determined, j is chosen as small as possible. I(x) i -ε(I i ; B)(x) is added to the component j. More formally, the erosion can be written as:

           ε i (I; B)(x) i = ε(I i , B)(x) ε i (I; B)(x) j = I j (x) + I i (x) -ε(I i , B)(x) for j = min(arg min k =i (δ(I k , B))) ε i (I; B)(x) k = I k (x) for k = i, j (3.19)
This rule is an interesting candidate of erosion ε i , however it does not provide idempotency of ϕ i and γ i . Indeed we did not nd a satisfying rule providing the stability of ϕ i and γ i without restricting the generality of the framework. The rst restriction consist in working only on the edges of the simplex, that is to say with distribution composed of two components at most, i.e., only two values of the mixture (a 1 , • • • , a n ) can be dierent from 0. If this set of distribution can not longer deal with triple points, it is still enough to describe transitions between two objects. The second restriction makes the lling rule depending on the location on the image. Consider the following rule. A grey-scale erosion of the component i is performed. When at x, i is the majority component and when there exists a component j = i majority at p ∈ B x , the value of the component j becomes 1 -ε(I i , B)(x) and all other components are set to 0. If i is not the majority component at x, or if there is no other majority component in B x , the remaining component are rescaled.

Formally, the restricted erosion is given by

                     ε i (I; B)(x) i = ε(I i , B)(x) if I i (x) > 0.5 and max k =i (δ(I k , B)) > 0.5 ε i (I; B)(x) j = 1 -ε(I i , B)(x) for j = min (arg max k =i (δ(I k , B))) ε i (I; B)(x) k =i,j = 0 else, ε i (I; B)(x) k =i = 1-ε i (I;B)(x) i 1-I(x) i .I(x) k (3.20)
As mentioned in Section 3.2.3 and illustrated Fig. 3.5, the opening operator γ i = δ i • ε i tends to be more interesting than the closing in the n-ary framework. The intuitive considerations about the openings and closings are conrmed by the fact that morphological properties are easier to recover for openings than for closings. It was the case for the n-ary lter introduced in Eq.(3.11), and it is now the case for the stability properties. Proposition 3.5.1 In the case of the n-ary continuous dilation of component i as in (3.18) and the restricted n-ary continuous erosion of component i, we have the following stability of the n-ary continuous opening:

γ 2 i (I; B) = γ 3 i (I; B).
Proof Let I ε i j = ε i (I; B) j . The grey-scale morphology gives

I γ i i = I γ 2 i i . Let us rst show that for all j, if I γ i j (x) > 0.5 then I γ i j (x) = I γ 2 i j (x).
The result is clear for j = i. For j = i, since I γ i i (x) is not the majority component of I γ i (x), the j-th component is simply rescaled after ε i . Thus I ε i •γ i j (x) ≥ I γ i j (x) > 0. Since δ i aects the components dierent from i only by a rescaling, and that

I δ i •ε i •γ i i (x) = I γ i i (x), we have that I δ i •ε i •γ i j (x) = I γ i j (x).
We will now show that for I γ i j (x) ≤ 0.5, I

γ 2 i j (x) = I γ 3 i j (x). I ε i •γ i j (x)
is determined by the majority components of I γ i (p) for p ∈ B x . Since the majority components are stable between I γ i and I γ 2 i , we have that

I ε i •γ i j (x) = I ε i •γ 2 i j (x). Finally since δ i only introduce rescaling on component dierent from i, I δ i •ε i •γ i j (x) = I δ i •ε i •γ 2 i j (x)
. Hence the result.

Experimental results

Let

Conclusions and perspectives

We proposed here an approach to extend mathematical morphology to images composed of n independent labels (or classes). The approach presents two key particularities: rst, the number of elementary operators is increased; second, the absence of an indeterminate class. This absence of indeterminate label transforms the notion of erosion the problem into a problem of gap lling.

In our framework, classical theory of morphological adjunction does not apply and consequently, exploring properties of composed operators requires specic considerations. Hence, we proved that some of the elementary properties of standard morphological operators are preserved, such as the idempotence of openings and closings per label. Despite its quasi experimental validity, the main lost property is the granulometric semigroup. Beyond the mathematical properties, one of the natural consequences of this n-ary framework is the denition of a new reconstruction operator. The main application of the proposed operators is the ltering of image objects according to their size/shape, the presented examples demonstrate the relevance of the n-ary operators. Fuzzy classication is a natural extension of the discrete n-ary framework to the continuous case, by working on the edges of the (n-1)-simplex. Even if theoretical properties are more dicult to recover in continuous cases, the result obtained in practice are promising.

Our future research will focus on the improvement of the theoretical framework for the continuous cases, including in particular a relevant denition of geodesic reconstruction, as well as on the study of more practical case-studies and applications.

Part II.

Probability density estimation on

Riemannian spaces

Abstract

This short chapter introduces theoretical notions and notations necessary for a good understanding of the second part of the thesis. We chose not to introduce the notion of manifold.

Indeed, expect spaces of rotations, every space present in the following study can by described by a unique chart. In other words, they are homeomorphic to R n . There is then no need to introduce an abstract denition of the tangent space based on equivalence classes.

The reader interested on a deeper insight into measure theory is referred to [START_REF] Tao | An Introduction to Measure Theory[END_REF],

for more details on Riemannian geometry to [START_REF] Berger | A Panoramic View of Riemannian Geometry[END_REF] and for a reference on analysis on symmetric spaces to [START_REF] Helgason | Geometric Analysis on Symmetric Spaces[END_REF].

Basics on measure theory and probability

Let Ω = ∪ {ω} be a set of outcomes. Let P(Ω) be the set of subsets of Ω and A be a subset of P(Ω). An element of P(A) is called an event. A represents the subset of P(Ω) where it is possible to quantify the probability. A probability measure is a function µ Ω that goes from A to R + which associates a probability to each event. Measure µ Ω and set A must verify several properties. First, it is imposed that:

Ω ∈ A, µ Ω (Ω) = 1.
Let A, B be in P(Ω) such as A ∩ B = ∅. If it is possible to quantify the probability of A and B, that is to say A, B ∈ A then it is possible to quantify the probability of the union:

µ Ω (A ∪ B) = µ Ω (A) + µ Ω (B).
As a consequence A∪B ∈ A. The property is extended to any family for disjoint countable sets (A i ) i∈I⊂N :

∪ i A i ∈ A, 56
4.1. Basics on measure theory and probability

µ Ω (∪ i A i ) = Σ i µ Ω (A i ).
The last property required of A is:

∀A ∈ A, A ∈ A,
where A is the complement of A. Under these assumptions, A is called a σ-algebra and µ a probability measure on Ω. In a more general way, a measure µ is a probability measure without the constraint µ(Ω) = 1. In a probability space, the measure is used to quantify the probability of an event. And in a more general way, the measure is used to formalize the notion of area and volume.

We can now dene the notion of random variable. Let V be a space with its own σalgebra. Let A Ω and A V be respectively the σ-algebras of Ω and V . A random variable X is a measurable function from Ω to V . A measurable function is a function such that:

A ∈ A V ⇒ X -1 (A) ∈ A Ω .
The measurable hypothesis implies that the variable X can be seen as a transport of the measure of probability µ from the space Ω to the space V . Indeed, the function µ X on A V dened by:

µ X (A) = µ Ω (X -1 (A)), is a probability measure on V .
The concept of density can now be introduced. Let E be a set and A E a σ-algebra of E. Let µ 1 and µ 2 be two measures on E. A measure is a way of measuring volumes.

Consequently, measures play a central role in the notion of integral: the integral of an indicator function of a set A is equal to its volume, that is to say to its measure. The integral of an ordinary function is dened as the limit of a weighted sum of integral of indicator functions. The notion of integral depends then on the choice of a reference measure. We note f dµ the integral of the function f with respect to the measure µ. We say that f = dµ 1 dµ 2 is the density of the measure µ 1 with respect to the measure µ 2 if we have

∀A ∈ A E , µ 1 (A) = A f dµ 2 .
Let µ V be a reference measure on V . For some specic random variable X, the measure µ X has a density f with respect to the reference measure of V . Then f is called the density of the random variable X. In a nite dimensional vector space there is a unique measure, up to a scaling factor, invariant under translations. The translation invariant measure that normalizes the unit hyper-cube is called the Lebesgue measure. In the case where V is a vector space the reference measure µ V is often the Lebesgue measure. Fig. 4.1 illustrates a change of reference measure.

Let E and F be sets equipped with σ-algebras, and f : E → F an application . Any measure µ dened on E can be transformed into a measure f * (µ) on F by:

f * (µ)(A) = µ(f -1 (A)). (4.1)
The smallest σ-algebra of E × F that contains all the products U × V, (U, V ) ∈ A E × A F is called the product σ-algebra of E × F . A product measure µ on E × F is a measure such that there exist two measures µ E and µ F respectively on E and F with:

∀(A, B) ∈ A E × A F , µ(A × B) = µ E (A)µ F (B).
L(γ) = b a (φ -1 • γ) (t), (φ -1 • γ) (t) (φ -1 •γ)(t) dt, (4.2) where •, • (φ -1 •γ)(t)
is the scalar product attached to the point (φ -1 • γ)(t). The notion of shortest path between two points induces a distance on M. The distance d(p, q) associated with the Riemannian metric is given by d(p, q) = inf γ {L(γ)},

(4.3)
where the inmum is taken over all the smooth paths from p to q. A path realizing this minimum is called a geodesic path. Geodesic paths can be seen as straight segments on M. Geodesics are paths which are locally shortest paths.

Given p ∈ M with φ(x) = p, the set of vectors u ∈ R n attached to x is noted T p M and is called the tangent space at p. Unlike the case of manifolds where the tangent space is dened as an equivalence class isomorphic to R n , we simply have here T p M = R n .

It can be shown that given a < 0, 0 < b, vector u ∈ R n , there is only one geodesic γ, [a, b] → M such that γ(0) = p with tangent vector (φ

-1 • γ) (γ -1 (p)) = u such that ||(φ -1 • γ) (t)|| (φ -1 •γ)(t) = ||u|| x for all t ∈]
a, b[. It can be then shown that there exists a unique geodesic such that γ(0) = p, (φ -1 • γ) (γ -1 (p)) = u and that its domain can not be extended. This geodesic is noted geod p (u). The exponential map, see Fig. Note that this exponential has generally no link with the classical exponential application. However, they happen to coincide in specic cases.

The radius of injectivity r inj p of exp p is the largest r such that exp p restricted to the centered ball of radius r is injective. Inside this ball, log p (.) denotes the inverse of exp p (.).

Let S be a subset of M homeomorphic to an open subset of R m with m ≤ n. The Riemannian metric on M naturally induces a Riemannian metric on S. The metric on S is called the induced metric. In the expression of the geodesic distance d S (p, q) = inf γ {L(γ)}, the inmum is then taken over all the paths γ that stay in S.

A Riemannian metric also induces a measure. The matrix of the metric G(x) is symmetric and can be written G(x) = AA t where A is a n × n inversible matrix and represents a local linear change of coordinates that induces the new scalar product G(x). The modication of a unitary volume associated with the linear change of coordinates is given by det(A) = det(G(x)).

(4.4)

The density of the Riemannian measure with respect to the Lebesgue measure of U φ is given by x ∈ U φ → det(G(x)).

In many situations, given a metric, it is impossible to obtain an analytic expression of the distance between two points of M . An isometry of M is bijection h : M → M such that d(p 1 , p 2 ) = d(h(p 1 ), h(p 2 )). Using the parametrization φ an isometry is indierently seen as an application h : M → M or h : U φ → U φ . The set ISO(M, G) of isometries is a group: i) the composition of two isometries is an isometry and since an isometry is a bijection, ii) every isometry has an inverse which is also an isometry. Many properties of the geometry induced by the metric can be deduced from the study of ISO(M, G). A group action from a group A on a set B is an application act : (A, B) → B compatible with the left-multiplication of A:

∀a 1 , a 2 ∈ A, b ∈ B, act(a 2 , act(a 1 , b)) = act(a 2 a 1 , b).
The orbit of an element b ∈ B is dened as:

Orb(b 0 ) = {b ∈ B |∃a ∈ A, act(a, b 0 ) = b}.
An action is said transitive when Orb(b) = B. In other words, from a point b in B, it is possible to reach every point of B through the action act. Note that if Orb(b) = B for one b in B, it is also true for all the points in B.

Group of isometries ISO(M, G) naturally acts on M: act(iso, p) = iso(p), iso ∈ ISO, p ∈ M.

Riemannian manifold M is said to be homogeneous if the action of ISO(M, G) is transitive, i.e., given two points of p 1 , p 2 ∈ M, there is always an isometry iso such that p 2 = iso(p 1 ).

In a homogeneous space, all the points are similar in term of geometry. Vector spaces and spheres are examples of homogeneous spaces. However an open disc is not a homogeneous space: if p and q are two points such that d(0, p) is dierent from d(0, q), there is no isometry iso disc of the disc such that q = iso disc (p).

Let p be a point in M. The group K p of isotropy of p is the subgroup of ISO(M, G) that xes p, i.e.,

K p = {k ∈ ISO(M, G)|k(p) = p}.
M is homogeneous if and only if M is in bijection with the quotient space ISO(M, G)/K p .

A space M is said to be symmetric if for all p ∈ M there exists an element k ∈ K p such that:

dk = -id,
where dk is the dierential of k. In other words, reversing the geodesics starting from a point p is an isometry: central symmetries are isometries. It can be shown symmetric spaces are homogeneous.

Convolution on Riemannian spaces

Space M is said to be isotropic if for all p ∈ M and for all couples of vectors u, v ∈ R n , there is an element k of the group of isotropy of p such that:

dk(u) = v.
At each point p, an isotropic space M looks the same in every direction. Very few spaces are isotropic. However, several spaces like torus are locally isotropic. It can be shown that isotropic spaces are symmetric and therefore homogeneous.

Convolution on Riemannian spaces

We introduce now some basic ingredients of analysis on Riemannian spaces. The notion of convolution f * g of two functions M → R is highly related to the notion of homogeneity. In R n the convolution is dened as follows:

(f * R n g)(p) = R n f (x)g(p -x)dx, p ∈ R n . (4.5)
Thus, the convolution of a function f by a kernel g consists in the integral of the translated kernel in each point x of the support space, weighted by f (x). The group law + of R n is an isometry that enables to transport the kernel in the whole space.

If the space is not homogeneous, the transport of the kernel might be impossible without modication of its shape.

Examples. Let g 0 be the triangular kernel centered in 0 dened on the open unit disc of R by g(x) = 1 -|x|, x ∈ [-1, 1]. Kernel g cannot be naturally transported at any other point of the disc. On the other hand if g 1 is dened on the unit circle of R 2 by g 1 (e iπθ ) = 1 -|θ|, θ ∈ [-1, 1] can easily transported to any point p = e iπθ 0 of the circle: g p (e iπθ ) = g 1 (e iπ(θ-θ 0 ) ).

In an homogeneous space, it is possible to transport a kernel from a reference point to any other point by an isometry. However, there are often several isometries that send a point p ref on a point q. If iso ∈ ISO(M, G) is such that iso(p ref ) = q, then for any k ∈ K p ref , iso(k(p ref )) = q. Since there are several potential candidates to transport the kernel from p ref to q, the convolution kernel g is required to be invariant under the isotropy group of the reference point. Then the transport of the kernel from p ref to q is independent of the choice of the isometry. Let iso q,p ref be any isometry that sends q on p ref such that g q (p) = g p ref (iso q,p ref (p)). Thus in an homogeneous space M, the convolution of a function f ∈ L 2 and a K p ref -invariant kernel g p ref ∈ L 2 is naturally dened by:

(f * M g)(p) = q∈M f (q)g q (p)dµ, (4.6)
where µ is the Riemannian measure. Note that the convolutions * M , and the standard convolution on R n slightly dier. Indeed, on R n the isometry used to send p on q is always the translation of vector q -p so that the classical convolution on R n does not require an isotropic kernel.

In R n the Laplacian of a function f is dened as follows:

∆ R n f = div(gradf ) = ∂ 2 i f.
It is possible to adapt the denition of the gradient and of the divergence operators on (M, G). The expression of the so-called Laplace-Beltrami operator becomes:

∆ M f (x) = 1 det(G(x)) i,j ∂ ∂x i ( det(G(x))g i,j ∂f ∂x j ),
where g i,j are the coecient of G(x) -1 . In R n , the functions involved in the Fourier transform are eigenfunctions of the Laplacian operator:

∆ R n e i ω,. = ||ω|| 2 e i ω,. , ω, x ∈ R n .

The generalization of the Fourier transform on symmetric spaces is called the Fourier-Helgason transform. As trigonometric exponentials, the Helgason functions e ω (x) are eigenfunctions of the Laplace-Beltrami operator. The theoretical construction of Helgason functions e ω (x) (generalization of plane waves) can be found in [START_REF] Helgason | Non-Euclidean Analysis[END_REF]. Explicit expression of these functions are provided when necessary in next chapters. The transform and the inverse transform are of the form:

H(f )(ω) = x f (x)e ω (x)dµ 1 (x), (4.7) f (x) = ω H(f )(ω)e ω (x

Image histogram and density estimation

In this section we set the theoretical links and notation between image histograms and density estimation. Let us consider an image I as the map:

I : Ω → V p → I(p)
We have for instance V = R for grey-scale images or V = R n for multispectral images. Domain Ω is the support space of pixels/voxels, typically a subset of R 2 or R 3 such as a rectangle or a parallelepiped.

We endow Ω with a measure µ Ω that reects its geometry: since the spatial distances are usually evaluated using Euclidean distance, this measure is typically the Lebesgue measure. Indeed the prior importance of an area is often proportional to its Lebesgue measure. Measure µ Ω is generally nite, that is to say µ Ω (Ω) < ∞.

The application I transports the measure µ Ω on V . The image measure on V is noted I * (µ Ω ), see Eq. (4.1). The measure I * (µ Ω ) represents the distribution of values. It often carries a useful synthetic information of the image I. Thus, one might be interested in the study of this distribution. We have access to a nite number of evaluations of I in points p ∈ Ω called pixels, from which we try to deduce information on I * (µ Ω ). The value space V is usually endowed with a reference measure µ V . The standard strategy is to build a density f with respect to µ V from the set of observations {I(p), p a pixel } such that the associated measure approximates I * (µ Ω ). Note that in many situations I * (µ Ω ) does not have a density with respect to µ V . This is for instance the case for color images I : [0, 1] 2 → R 3 when µ Ω and µ V are Lebesgue measures. Then the density f does not represent exactly the measure I * (µ Ω ).

We assume that pixels p are uniformly distributed with respect to µ Ω . Roughly, for A ⊂ Ω,

1 N p 1 A (p) ≈ µ Ω (A),
where N is the number of pixel and 1 A the indicator function of A. Then, if µ Ω is normalized, the set of pixels can be considered as a set of N independent realizations of a random variable of law µ Ω . It follows directly that the set {I(p), p a pixel } can be considered as a set of N independent realizations of a random variable of law I * (µ Ω ). The study of I * (µ Ω ) is thus often interpreted as a problem of probability density estimation. This is for instance the case of color histograms developed in Section 5.3. However, note that this way of addressing the problem do not take advantage of the specicity of the situation.

Let X be a random variable from (Ω, µ Ω ) to (V, µ V ). Assume known a set of independent realizations X(ω i ), {ω i ∈ Ω}. The estimation of the density of X * (µ Ω ) with respect to µ V is a very dicult problem since no assumption is made on X(ω / ∈ {ω i }). The amount of information contained in {X(ω i )} about X * (µ Ω ) is extremely weak, on might even say null. Thus the estimated density cannot be determined with condence, but can only be guessed. When building an histogram, the choice of the bin size is an essential parameter and there is no precise way of determining if two peaks should be gathered in the same bin.

However, the situation in image processing is dierent. Unlike the standard case of random variables, the space Ω often has a metric structure and the function I has regularity properties with respect to this metric. Thus, knowing the values {I(p i )} brings information on I(p / ∈ {p i }). The set {I(p i )} carries then a signicant amount of information about I * (µ Ω ). Assume for instance a rule of interpolation of I(p / ∈ {p i }) from {I(p i )}. The density f of I * (µ Ω ) with respect to µ V can then be exactly computed using formulas such as the coarea formula, see [START_REF] Morgan | Geometric measure theory: a beginner's guide[END_REF].

We rst give an almost rigorous formula for I : [0, 1] 2 → R, where spaces are endowed with Lebesgue measures. When possible let γ h (t) be a unit speed parametrization of

I -1 (h). Then we have f (h) = dI * (Leb [0,1] 2 ) dLeb R (h) = 1 ||∇I(γ h (t))|| 2 dt.
For a function φ, R n → R k , the k-dimensional Jacobian determinant J k (x) is dened by the maximal volume of the image of a unit k-dimensional cube by the dierential of φ. The general formula for I,

[0, 1] n → R k with k < n is, f (h) = dI * (Leb [0,1] n ) dLeb R k (h) = I -1 (h) 1 J k (x) dH n-k (I -1 (h)),
where H n-k refers to the (n -k)-Hausdor measure, see [START_REF] Morgan | Geometric measure theory: a beginner's guide[END_REF].

Note that an interpolation of I : [0, 1] 2 → R by pieces of linear functions leads to an elementary computation of the density estimator f . The previous problem of the size of bins no longer exists, the shape of the computed density f is intrinsic to the interpolated image. While they did not follow exactly the same line of thoughts, Authors of [START_REF] Duy | Integrating Isosurface Statistics and Histograms[END_REF] ended up at the same conclusions. Their work still remains quite undisclosed in the image processing community. Fig. 5.1 presents histograms of a 16 bits gray-scale image computed using dierent methods.

Perceptual color histograms

One sometimes considers the color space C as a part of a three dimensional Euclidean space R 3 . Under this assumption, the histogram of a color image can be built in the same way as for gray-scale images. However, the distances induced on colors by the human perceptual system cannot be represented by a Euclidean space structure. Observation showed that the perceptual relation between colors is better represented in the framework of Riemannian manifolds.

The local metrics of the Riemannian structure are experimentally measured by a set of ellipses, such as the MacAdam ellipses introduced in Macadam (1942), BFD-P by [START_REF] Luo | Chromaticity-discrimination ellipses for surface colours[END_REF] and RIT-DuPont by [START_REF] Berns | Visual determination of suprathreshold color-dierence tolerances using probit analysis[END_REF]. This Riemannian structure makes the construction of the histogram dicult. On the one hand, except rare situations, there are no regular tilings of the space. On the other hand, kernel methods have been generalized to Riemannian manifolds in [START_REF] Pelletier | Kernel density estimation on Riemannian manifolds[END_REF], but requires the knowledge of the geodesics.

We propose here an approach that takes into account the Riemannian structure of C while keeping the computation in the Euclidean framework. Thus we propose a way of building histograms for images I : [0, 1] 2 → C that respects the perceptual distances better than histograms built in Euclidean spaces, and without increasing the computation time.

Perceptual metric on colors

Already Riemann used colour as an illustration of the applicability of his geometry in Riemann (1867), and concrete examples of such colour geometries were developed by Helmholtz (1891), [START_REF] Schrödinger | Grundlinien einer Theorie der Farbenmetrik im Tagessehen[END_REF] and [START_REF] Stiles | A modied Helmholtz line-element in brightness-colour space[END_REF]. A scalar product is entirely determined by its unit ball. Expressed in vector coordinates, the associated unit ball takes the form of an ellipse in two dimensions or of an ellipsoid in three dimensions. Thus, the Riemannian metric is given by a eld of ellipses or ellipsoids.

Ellipses, local metric

The rst experimental determination of the eld of ellipses describing the Riemannian metric of the colour space was performed by [START_REF] Macadam | Visual Sensitivities to Color Dierences in Daylight[END_REF]. The experiment consisted of about 25 000 colour matches with one observer, and the ellipses were derived from the covariance matrices of the repeated observations. Later, it has become common practice to denote ellipses obtained in this manner as JND (just noticeable dierence) ellipses or ellipsoids.

Later, another type of experiment has become more commonplace. Pairs of colours that are barely perceptually dierent, are presented to the observer, who is given the task to estimate the magnitude of the perceptual distance using a set of standard pairs. Ellipses, ellipsoids and metrics obtained in this way are normally denoted supra-threshold ellipses.

Examples of supra-threshold color dierence based data include BFD-P, see [START_REF] Luo | Chromaticity-discrimination ellipses for surface colours[END_REF] and RIT-DuPont, see [START_REF] Berns | Visual determination of suprathreshold color-dierence tolerances using probit analysis[END_REF].

Global model

Data sets of measurements provide information on distances through the local metric or through distances between specic colors. A global model provides an analytic expression of the distance between two arbitrary colors. The closest the proposed expression is to the Riemannian perceptual distances, the better the model is. The more conventional procedure for going from a tristimulus space to a space closer linked to a perceptual homogeneous space typically includes the following steps.

First, apply a linear transform in the tristimulus space such that the base gets close to the cone fundamentals of the retina.

Secondly, perform a non-linear compression of the coordinates (e.g., logarithmic or cubic root) in order to mimic the non-linear response of the human visual system.

Finally, perform a linear transformation of the resulting coordinates in order to correspond better to the perceptual attributes of color.

Typically, the rst coordinate is a weighted sum of the coordinates and represent a lightness correlate, whereas the two other coordinates are weighted dierences, and represent color opponent channels such as, e.g., redgreen and blueyellow.

In order to identify the dierent parameters of the various transforms, dierent optimisation criteria are used. In the CIELAB colour space, see [START_REF] Robertson | The CIE 1976 Color-Dierence Formulae[END_REF], the parameters were optimised so that the lightness should correspond to perceived lightness, and that the Euclidean metric in the resulting space should correspond to perceptual colour dierences.

For the IPT colour space, see F Ebner (1998), the parameters were optimised in order to achieve a constant perceived hue along straight radial lines in cylindrical coordinates. It is furthermore reasonably well established that in such perceptual spaces, the Euclidean metric is not the one best corresponding to the perceived colour dierences, and other models have been proposed, see, e.g., [START_REF] Luo | The development of the CIE 2000 colour-dierence formula: CIEDE2000[END_REF] and [START_REF] Farup | Hyperbolic geometry for colour metrics[END_REF]. In the hyperbolic models proposed in [START_REF] Farup | Hyperbolic geometry for colour metrics[END_REF], histograms can be computed using adapted kernels, see next Chapter.

Kernel density estimation

There are various ways of addressing the problem of probability density estimation. In the Euclidean context the most popular techniques are mainly the histograms, the kernels, and the characteristic function density estimator. The characteristic function density estimator consists in the estimation of the Fourier transform or series of the density. In what follows, we chose to focus on the kernel method. Recall that the kernel method in the Euclidean case has the following expression for the density estimator:

f (x) = 1 k p i ∈{pixels} 1 λ n K ||x -I(p i )|| λ , (5.1)
where λ is a scaling parameter, n the dimension of the space, k the number of pixels, and K : R + → R + a map which obeys the following properties:

-

R n K(||x||)dx = 1, - R n xK(||x||)dx = 0, -sup(K(x)) = K(0).
In the present work, we assume a supplementary condition of bounded support K(x > 1) = 0.

Local Euclidean approximations

In general Riemannian manifolds, computing the distance between two arbitrary points given the metric eld is a dicult problem. Indeed, nding the distance is a minimization problem over a set of paths. However, for two close points, the local metric provides a satisfying approximation of the Riemannian distance. A probability density measures the ratio between the probability of an innitesimal volume element and its volume. It is thus a local notion. The central idea of our contribution is to take advantage of the fact that histograms mainly involves local phenomena. Since in a Riemannian manifold the computation of an histogram does not involve computation of long geodesics, the need of a global model that provides distances between every pairs of colors is of lower importance than in most applications. Let us assume that these ellipses represent the local perceptual metric. Let c be a point where the metric has been measured through the ellipse E c . In a neighborhood of c, computing distances using the metric measured at c is a better approximation of the perceptual distance than using the canonical Euclidean distance of the ab plane. At a point p where the metric is originally unknown, a metric interpolated from the neighbor points c i has all the odds of being more relevant than the canonical Euclidean metric of the map, see Fig. 5.2(c).

Let d R (p, q) be the perceptual distance between color p and color q. d R (p, q) is the Riemannian distance induced by the eld of ellipses. Let ||p-q|| be the distance associated with the canonical scalar product of the ab plane, and ||p -q|| c be the distance associated with the scalar product induced by the ellipse E c . Let B(c, R) and B c (c, R) be the respective balls of center c and radius R. The previous discussion can be formalized as follows. It can be shown that:

lim x→c ||x -c|| c d R (x, c) = 1, while if ||.|| c = ||.
||, the equality case being exceptional, thus

lim x→c ||x -c|| d R (x, c) = 1.
Therefore for such a c there exists A > 0 such that,

∀R > 0, ∃x ∈ B(c, R), A < ||x -c|| d R (x, c) -1 .
(5.2)

On the other hand there exists a real positive number

R c = R c,A such that, ∀x ∈ B(c, R c ), ||x -c|| c d R (x, c) -1 < A.
(5.3)

We have

sup B(c,Rc) ||x -c|| c d R (x, c) -1 < A < sup B(c,Rc) ||x -c|| d R (x, c) -1 . Thus for x ∈ B(c, R c ), ||x -c|| c is preferred to ||x -c||. Consider a kernel K and a scaling parameter λ such that λ ≤ R c and B c (c, λ) ⊂ B(c, R c ). For x ∈ B(c, R c ), K ||x-c||c λ is preferred to K ||x-c|| λ . For x / ∈ B(c, R c ), K ||x-c||c λ = K ||x-c|| λ = 0.
Therefore, under these assumptions on the scaling parameter λ, the histogram estimator

f (x) = 1 k p i ∈{pixels} 1 λ n K ||x -I(p i )|| I(p i ) λ ,
(5.4) is preferred to the classical histogram. We think that the hypothesis on λ is reasonable in practice, its validation is a subject of further research. Note that the higher the resolution of the image is, the smaller λ is and then the more the hypothesis becomes reasonable.

Metric interpolation and Euclidean approximation

Let M be topological space, homeomorphic to an open subset of R n . Let φ be an homeomorphism from U φ ⊂ R n to M. A set of scalar products G c i is given for a set of points {c i } ∈ M. We consider here the problem of interpolation of the eld of metrics. Let F 1 and F 2 be two smooth metric elds that coincide with the observed ellipses at the points {c i }. Despite the intuition, if no assumption is made on φ regarding the Riemannian distance, there are no criteria that enables to prefer F 1 or F 2 . The problem of metric tensor interpolation is thus a dicult problem. In this work, we adopt an elementary solution.

Ellipses are represented in the projective ab plane. A Delaunay triangulation with respect to the canonical Euclidean metric of the plane is performed on the set {c i }, Fig. 5.2(b). At a point p in the triangle c i c j c k the parameters of the interpolated ellipse E p are linearly interpolated between the parameters of E c i , E c j , and E c k with respect to the barycentric coordinates of p. If p does not belong to one of the triangles of the Delaunay triangulation, we set E p = E q where q is the projection of p on the convex hull of the set of centers.

Experimental results

The RIT-DuPont dataset from [START_REF] Berns | Visual determination of suprathreshold color-dierence tolerances using probit analysis[END_REF] shows that the perceptual metric is Let us consider the color image given in Fig. 6.2.2(a) and the three alternative color histograms compared in Fig. 6.2.2. Histograms (b) and (c) aim a studying the density f with respect to the perceptual Riemannian volume measure. The main dierence between (c) and (d) is that in (d) the shape of the kernel follows the Riemannian metric. The density with respect to the Euclidean measure is visibly dierent from the histogram with respect to the Riemannian measure. The amplitude of the upper spot, representing white colors, is signicantly decreased when using the Riemannian measure. Perceptually, this results from the fact that the eyes have an higher sensitivity around white than around blue.

Conclusion and perspectives on perceptual color histograms

Given a set of ellipses representing the perceptual metric on colors, we proposed an approach for histogram computation that takes into account the Riemannian structure of the perceptual metric without introducing supplementary computational complexity. Indeed, the step of ellipses interpolation only has to be achieved once and does thus not introduce computational complexity. The relevance of the approach is conditioned by the relevance of the set of perceptual ellipses and the quality of the interpolation. The deep problem of metric tensor interpolation has been partially left aside and will be subject of future research. The second topic of our future research will be on the convergence of the proposed histogram to the density of the underlying random variable with respect to the interpolated Riemannian metric. The most common geometries of Gaussian laws are the Fisher information metric and the Wasserstein metric. The Fisher metric on multivariate centered Gaussians happens to be similar to the ane invariant metric on SPD(n), that is to say the metric invariant under the action of invertible matrices Gl(n), i.e.,

G • Σ → GΣG t .
Furthermore the Fisher metric on univariate Gaussians is a hyperbolic metric of dimension two. Thus understanding the Fisher metric has consequences beyond Gaussian laws. The Wasserstein metric is another name for the earth mover's distance. Let S be a set of probability distributions on E. A fundamental dierence between the Fisher metric and the Wasserstein metric on S, is that the Fisher metric do not relay on a geometry of E while the Wasserstein metric does. Indeed, the Wasserstein metric depends on the geometry of E through a transportation cost. The dependence on the geometry of E constitutes one of the main specicities and interests of the Wasserstein metric.

Median and means are fundamental quantities in signal and image processing. Several studies already address their denition and computation for the Fisher and the Wasserstein metrics, see for instance [START_REF] Pennec | Probabilities and Statistics on Riemannian Manifolds : A Geometric approach[END_REF]; [START_REF] Yang | Riemannian Median and Its Estimation[END_REF]; [START_REF] Agueh | Barycenters in the Wasserstein Space[END_REF][START_REF] Arnaudon | Riemannian Medians and Means With Applications to Radar Signal Processing[END_REF]; [START_REF] Bini | Computing the Karcher mean of symmetric positive denite matrices[END_REF]. Although being a secondary problem, density estimation is an important tool of signal and image processing which is still little studied on spaces of Gaussian laws. It is in particular useful, for instance, for segmenting point clouds or for Bayesian classication. The study of distributions on SPD(n) is subject to recent studies. In [START_REF] Said | Riemannian Gaussian Distributions on the Space of Symmetric Positive Denite Matrices[END_REF], it is proposed for instance a generalization of the Gaussian distribution on SPD(n). To our knowledge, the existing literature dedicated to density estimation is mainly restricted to Huckemann et al. ( 2010), Peter T. [START_REF] Kim | Deconvolution Density Estimation on the Space of Positive Denite Symmetric Matrices[END_REF] and [START_REF] Asta | Kernel Density Estimation on Symmetric Spaces[END_REF]. These previous works have focussed on techniques derived from the orthogonal series density estimation for the Fisher metric in the case of centered multivariate and univariate Gaussian laws.

This chapter is a practical guide to density estimation on spaces on Gaussian laws and can be summarized as follows.

Main techniques of density estimation on Riemannian manifolds are reviewed in Section 6.2.

Application to each specic situation of Gaussian law spaces are discussed in Sections 6.4-6.13.

The computation of kernels for the Fisher and Wasserstein metric are provided. The expression of the dierent factors constituting these kernels in the case of the Fisher metric are already known in the literature, see [START_REF] Terras | Harmonic Analysis on Symmetric Spaces and Applications II[END_REF]; [START_REF] Wijsman | Invariant measures on groups and their use in statistics[END_REF].

The contributions of Sections 6.6.2 and 6.7.2 consist in presenting a computation relying on a low level Riemannian geometry and Lie groups, and providing a readyto-use expression of the kernels. To our knowledge, the expression of kernels for the Wasserstein metric is new.

Being an important problem for practical applications, the density estimation of partial quantities such as mean eigenvalues and rotation is discussed in Section 6.13.

Non parametric density estimation on Riemannian spaces

Let Ω be a space, endowed with a σ-algebra and a probability measure p. Let X be a random variable Ω → M with M homeomorphic to U an open subset of R n . The homeomorphism is noted φ. M is equipped with a Riemannian metric G. The Riemannian measure is called vol and the measure on M induced by X is noted µ X . We assume that µ X has a density, noted f , with respect to vol, and that the support of X is a compact set noted Supp. Let (x 1 , .., x k ) ∈ M k be a set of draws of X.

The Dirac measure in point a is dened as:

δ a (U ) : 1 if a ∈ U 0 if a ∈ U Let µ k = 1 k i δ x i
denotes the empirical measure of the set of draws. This section presents the main techniques of estimation of f from the set of draws (x 1 , .., x k ). The estimated density in x in M is noted fk (x) = f (x, x 1 , ..., x k ). Observe that fk (x) can be seen as a random variable. The relevance of density estimation technique depends on several aspects. When the space allows it, the estimation technique should not privilege specic directions or locations. This results in an isotropy and a homogeneity condition.

In the kernel method for instance, a kernel function K x i is placed at each observation x i .

Firstly, in order to treat directions equally, the function K x i should be invariant by the isotropy group of x i . Secondly, for an other observation x j , functions K x i and K x j should be similar up to the isometries that send x i on x j . These considerations strongly depend on the geometry of the space: if the space is not homogeneous and the isotropy group is empty, these indierence principles have no meaning. The convergence of the dierent estimation techniques is widely studied. Results were rst obtained in the Euclidean case, and are gradually extended to the probability densities on manifold, see [START_REF] Hendriks | Nonparametric Estimation of a Probability Density on a Riemannian Manifold Using Fourier Expansions[END_REF][START_REF] Pelletier | Kernel density estimation on Riemannian manifolds[END_REF][START_REF] Huckemann | Möbius deconvolution on the hyperbolic plane with application to impedance density estimation[END_REF]; [START_REF] Asta | Kernel Density Estimation on Symmetric Spaces[END_REF]. The last relevant aspect, is computational.

Each estimation technique has its own computational framework, which presents pro and cons given the dierent applications. For instance, the estimation by orthogonal series presents an initial pre-processing, but provides a fast evaluation of the estimated density in compact manifolds. If we have an estimator f Eucl of the density with respect to the Lebesgue measure of U , noted Leb U , it is possible to obtain an estimator of the density with respect to the Riemannian measure vol:

Estimation from the Euclidean structure of the parametrization

f = f Eucl dLeb U dvol = 1 det(G(x)) f Eucl . (6.1)
The estimation f is a probability density with respect to the Riemannian measure vol. However the estimation does not respect homogeneity and isotropy considerations. If f Eucl is constructed from an Euclidean kernel K, the Riemannian shape of K will dier from x i to x j in an uncontrolled way.

Histograms

We provide here an L 2 convergence rate of the histogram when the random variable X is valued on a Riemannian space, as formulated in Section 4.2, and whose support is compact.

We chose not to introduce the notion of manifold in the thesis. The result is however valid on Riemannian manifolds. Even if this results might already be known, we have not been able to nd it in the literature. The closest result we are aware of in can be found in Geroy (1974), where the author provides a condition on the uniform convergence of the histogram when X is valued in a metric space, without assumptions on its support.

Let (A k j ) 0≤j≤J k be a partition of U . Let α j be the number of draws that lays in A k j , i.e.,

α j = Card{x i ∈ A k j }. The histogram of the draws (x i ) is the function fk = 1 k j 1 vol(A k j ) α j 1 A k j , (6.2)
where 1 E is the indicator function of E. fk is an estimator of the density f . We suppose that:

i. The density f is α-Lipschitz with respect to the Riemannian metric;

ii. There exists a constant C such that, for all k and all i,j, one has

vol(A k j ) vol(A k i )
≤ C, and C will be called the homogeneity factor; iii. β k = max j (R k j ) goes to zero when k goes to innity, where 2R k j is the diameter of A k j .

Estimator fk (x) is seen as a random variable of the set of draws. The quantity

e = E Supp (f (x) -fk (x)) 2 dx = Supp E((f (x) -fk (x)) 2 )dx,
is the mean square error of the estimation. It can be developed as:

e = Supp (f (x) -E( fk (x)) 2 dx + Supp E((E( fk (x)) -fk (x)) 2 ) = Supp (f (x) -E( fk (x))) 2 + Supp V ar( fk (x))dx = e 1 + e 2 . For x in A k j the random variable kvol(A k j ) fk (x) = α j has a binomial distribution B(k, µ X (A k j )), hence E( fk (x)) = µ X (A k j ) vol(A k j ) , V ar( fk (x)) = kµ X (A k j )(1 -µ X (A k j )) (kvol(A k j )) 2 ≤ µ X (A k j ) kvol(A k j ) 2 . It follows that e 2 ≤ j A k j µ X (A k j ) kvol(A k j ) 2 dx = j vol(A k j ) µ X (A k j ) kvol(A k j ) 2 .
Call N k the number of bins and V the measure of the support vol(Supp). By the assumption

(ii), CN k vol(A k j ) ≥ V, therefore e 2 ≤ j µ X (A k j ) k V CN k = CN k kV .
Note that j µ X (A k j ) = 1 since µ X is a probability. Assuming that the A k j are connected for each j, there exists y k j in A k j such that

E( fk (x)) = µ X (A k j ) vol(A k j ) = 1 vol(A k j ) A k j f (x)dx = f (y k j ).
We have

e 1 = j A k j (f (x) -E( fk (x))) 2 dx = j A k j (f (x) -f (y k j ) 2 dx,
and since f is α-Lipschitz we obtain

e 1 ≤ j (α.R k j ) 2 vol(A k j ) ≤ (αβ k ) 2 V.
Finally we have

e ≤ CN k kV + (αβ k ) 2 V.
Assume now that there exist positive real numbers C 1 ,C 2 such that we have

∀j, C 1 (R k j ) n ≤ vol(A k j ) ≤ C 2 (R k j ) n .
When the radius R goes to zero, the volume of a ball B R becomes proportional to R n , see [START_REF] Gray | The volume of a small geodesic ball of a Riemannian manifold[END_REF]. Thus the existence of C 1 ,C 2 is an assumption on the isotropy of the partition. It follows that for all i, j

C 1 (R k j ) n C 2 (R k i ) n ≤ C, which implies that min j (R k j ) ≥ C 1 CC 2 1/n β k = γ 1/n β k . Hence N k C 1 γβ n k ≤ V, and e ≤ C kV × V C 1 γβ n k + (αβ k ) 2 V = A kβ n k + Bβ 2 k = F (β k ),
where A and B are two constants depending on the geometry and on the density f . Let β k be the value that minimize the right hand side of the above inequality. We have:

β k = n.A 2Bk 1 n+2
, and therefore the minimum of F is:

min(F ) = D.k 1 n+1 .
The obtained convergence rates are similar to the rates of Euclidean regular multidimensional histograms. The constant D depends on the homogeneity and the isotropy of the partition through C, C 1 and C 2 . In practice, nding a partition with satisfying factors of homogeneity and isotropy can be dicult. It can be noted that since the bins are xed, even with ideal factors, C = 1, C 1 = C 2 , the estimation can not be exactly invariant to isometries.

Compared to other methods, the main advantage of histograms is often the low algorithmic complexity and the simplicity of use. The dierent costs that have to be taken into account are: the construction of the partition and the cost c(x) of the search of the bin A k j such that x ∈ A k j . Function c(x) depends on the geometry of the space. The second technical aspect that has to be taken into account is the nature of the representation of the partition: list, graph, matrix, etc. Indeed for further processing, the histogram is often viewed as an image, its structure plays then an important role.

Kernel

This part is mainly based on the work by [START_REF] Pelletier | Kernel density estimation on Riemannian manifolds[END_REF]. Let K : R + → R + be a map which veries the following properties: i)

R d K(||x||)dx = 1, ii) R d xK(||x||)dx = 0, iii) K(||x|| > 1) = 0, iv) sup(K(x)) = K(0).
Let p ∈ M and x ∈ U with φ(x) = p. In the general case, given a point p ∈ M, exp p denes an injective application only on a neighborhood of 0 ∈ T p M. However in every studied situations of the present work, the parametrization induced by exp p is injective on the whole space. Recall that when T p M is endowed with the scalar product G(x), Euclidean distances to p corresponds to Riemannian distance to p. The Lebesgue measure of T p M is noted Leb p . The function θ p dened by: θ p : q → θ p (q) = dvol dexp * (Leb p ) (q), Lemma 6.2.3 If φ is an isometry of M, then θ p (q) = θ φ(p) (φ(q)).

Recall that µ k is the empirical measure of the set of draws.

Theorem 6.2.4 If M is homogeneous, then fk = µ k * K.

Proof Since d(x i , .) and θ p (.) are invariant by the isotropic group of p ref , also is K. Thus the convolution is well dened. Using notations of Chapter 4, one has

(δ i * K)(p) = δ i (q) Kq (p)dµ = Kx i (p) = K(iso x i ,p ref (p)) = 1 kr n (θ p ref (iso x i ,p ref (p))) -1 K(d(p ref , iso x i ,p ref (p))) = 1 kr n 1 θ x i (p) K(d(x i , p)).
Hence,

µ k * K = i δ i * K = 1 kr n i (θ x i (.)) -1 K(d(x i , .)) = fk .
Therefore, when M is homogeneous, the density estimation proposed by [START_REF] Pelletier | Kernel density estimation on Riemannian manifolds[END_REF] is a convolution. Each location and direction are processed as similarly as possible. As a result, it can be easily shown that for an isometry α we have: fk (x, x 1 , ..., x k ) = fk (α(x), α(x 1 ), ..., α(x k )).

It must be noted that the kernel method requires the existence of explicit expression of distances and of the function θ p , which renders it unusable in several situations. The kernel density estimation presents better convergence rates than histograms. The only computational aspect that has to be taken into account is the cost of evaluation of the density at a point x in M. The computational cost of fk (x) is in O(k). The constant mainly depends on the computational cost of distances. Note that when the number of evaluations is important, histograms present lower computational cost.

The k-nearest neighbor is an interesting variant of the standard kernel density estimation. It consists in making the parameter r depend on x i , by setting r as the distance to the k-th nearest neighbor of x i . This estimator has been studied on Riemannian manifolds in [START_REF] Henry | k-Nearest neighbor density estimation on Riemannian Manifolds[END_REF].

Orthogonal series

Instead of directly studying the density f , the estimation is made from the estimation of the scalar product between f and a set of orthonormal functions {e j }. Recall that for g and h in L 2 (M):

g, h = M ghdvol,
where x denotes the complex conjugate of x. We have that:

f, e i = f e j dvol = E (e j (X)) .

Thus f , e j can be approximated by an estimation of the expectation, i.e., E (e j (X)) ≈

1 k k j=1 e j (x i )
Or in other words: e j f dvol ≈ e j dµ k . Now, given an integer N , the orthogonal series estimator is dened as:

f = 1 k N j=-N k i=1 e j (x i ) e j . (6.5)
The formula becomes an integral when the number of functions e j is uncountable. First the base has to be ordered such that the norm of the rest of the decomposition, i.e., |j|>N f, e j e j , decreases as fast as possible for regular functions. Second, in order to process locations and directions indierently, the basis functions must present regularity properties regarding the metric. Eigenfunctions of the Laplacian operator are very good candidates. Indeed the Fourier transform on R n , or the FourierHelgason transform on symmetric spaces, are highly related to the convolution by isotropic kernel. When the e j functions are Fourier functions, the estimation technique is also called the characteristic function method.

When the space is compact, the spectrum of the Laplacian operator is discrete. Under the compact assumption, let {e j } be an orthonormal Hilbert basis given by the eigenfunctions of the Laplacian. Given T > 0 we dene:

F T = {e j |∆e j = λe j with λ < T }.
Using this orthogonal series, the density estimator becomes:

f = 1 k e j ∈F T k i=1
e j (x i ) e j . [START_REF] Hendriks | Nonparametric Estimation of a Probability Density on a Riemannian Manifold Using Fourier Expansions[END_REF] provides the following convergence rates for density estimators based on eigenfunctions of the Laplacian.

Classical metrics on Gaussian laws

The hypothesis on σ is slightly dierent but highly related to the dierentiability hypothesis of the kernel method. Similar rates on the hyperbolic space of dimension two can be found here [START_REF] Huckemann | Möbius deconvolution on the hyperbolic plane with application to impedance density estimation[END_REF].

The FourierHelgason method on SPD(n) and on the hyperbolic space are convolutions:

fk = µ k * H -1 (1 F T ).
Recall that the FourierHelgason of the convolution of a function by a kernel is a product of transforms when the kernel is invariant by the isotropy group. This can be veried on expressions in following sections, by showing that H -1 (1 F T ) is invariant by the isotropy group SO(n) for SPD(n) and SO(2) for the hyperbolic space. Thus, as for the kernel method, each location and direction are processed as similarly as possible. The main drawback of this method in the non-compact case is its computational cost: the estimation requires the computation of an integral.

Curse of dimensionality

Observe that the convergence rates deteriorate with the dimension in k a/(b+n) where a and b are constants and n is the dimension of the space. To reach the same precision, the number of draws (x i ) has to increase exponentially with the dimension. Recall that the dimension of the space of centered n-dimensional Gaussian laws is of dimension 1 2 (n 2 + n) and that the dimension of n-dimensional Gaussian laws is of dimension 1 2 (n 2 +3n). Thus the quality of the density estimation quickly deteriorates with the dimension of the underlying space.

Fortunately, as shows section. 6.4.1, classical metrics on spaces of Gaussian laws allow to reduce the dimension by separating the dierent components of the Gaussian law (mean, eigenvalues and rotation). Dimension can also be reduced using principal geodesic analysis on manifolds, see [START_REF] Fletcher | Principal Geodesic Analysis on Symmetric Spaces: Statistics of Diusion Tensors[END_REF].

6.3. Classical metrics on Gaussian laws 6.3.1. The Fisher metric

We consider here a set M of probability measures on a measure space (X, σ X , µ) homeomorphic to an open subset of R n . Let θ = (θ 1 , ..., θ n ) ∈ U ⊂ R n be an homeomorphic parametrization of M. The distributions are assumed to have a density f (., θ). The Fisher matrix G(θ) associated to the parametrization θ is dened as follows, see Amari et al.

(

): g(θ) i,j = E ∂ln(f (x, θ)) ∂θ i ∂ln(f (x, θ)) ∂θ j = x ∂ln(f (x, θ)) ∂θ i ∂ln(f (x, θ)) ∂θ j f (x, θ)dx. (6.6) 2007
The matrix G(θ) is symmetric positive denite. The matrices {G(θ), θ ∈ U } induce a smooth eld of positive denite matrices on U , that is to say a Riemannian metric. It can be shown that the distances and the geodesic paths induced on M by this Riemannian metric do not depend on the choice of the parametrization. It is interesting to note that the Fisher metric do not rely on a metric on X.

Let P and Q be two probability measures. The relative entropy is dened as follows:

D(P, Q) = x ln dP dQ (x) dP.
Function D(., .) is not symmetric and is therefore not a distance. The relative entropy is also called the KullbackLeibler divergence. The KullbackLeibler divergence measures an informational dierence between P and Q.

Let D(θ, θ 0 ) be the relative entropy between the laws parametrized by θ and θ 0 . It can be shown that:

∂D(θ 0 , θ) ∂θ i (θ 0 ) = 0,
and

∂ 2 D(θ 0 , θ) ∂θ i ∂θ j (θ 0 ) = g(θ 0 ) i,j .
Thus, the Fisher matrix is also the Hessian of the KullbackLeibler divergence.

The Wasserstein metric/ Earth mover's distance

The Wasserstein metric on probability distributions is the metric of optimal transport as founded by [START_REF] Monge | Mémoire sur la théorie des déblais et des remblais[END_REF], see [START_REF] Villani | Optimal Transport[END_REF] for a modern and complete mathematical overview of the problem. If the probability distributions are seen as earth heaps, the Wasserstein distance between two distributions is the minimum cost that is needed to transform the rst heap into the second one. The cost of transport of earth depends of the amount of earth and the traveled distance. Due to this interpretation, the Wasserstein metric is also called the earth mover's distance.

More formally, let X and Y be two random variables valued in a space S endowed with a σ-algebra A S . Let µ X and µ Y be the measures induced on S by X and Y . Let Γ(µ X , µ Y ) be the set of measures on S ×S such that µ X and µ Y are the rst and the second marginals, i.e., for γ ∈ Γ(µ X , µ Y ), one has γ(A, S) = µ X (A), γ(S, B) = µ X (B), for all A, B ∈ A S . Thus Γ(µ X , µ Y ) represents the set of possible transports of the measure µ X on the measure µ Y . The p-Wasserstein distance between X and Y is then dened by:

d(X, Y ) = inf γ∈Γ(µ X ,µ Y ) S×S d(a, b) p dγ 1/p . (6.7)
For most families probability distributions, the Wasserstein distance has no explicit expression. Fortunately, the expression of the 2-Wasserstein distance between Gaussian laws has a simple form. Given two Gaussian distributions N 1 , N 2 of mean m 1 , m 2 and covariances Σ 1 , Σ 2 the distance becomes:

d(N 1 , N 2 ) 2 = |m 1 -m 2 | 2 + tr(Σ 1 ) + tr(Σ 2 ) -2tr( Σ 1/2 2 Σ 1 Σ 1/2 2 ).
Takatsu (2011) has showed that the 2-Wasserstein distance on Gaussian laws is a Riemannian metric distance. We restrict the present study to the case of the 2-Wasserstein distance.

The space of multivariate Gaussian laws

In this section we study the space of multivariate Gaussian laws under two classical Riemannian metrics: the Fisher metric induced by the Fisher matrix, and the Wasserstein metric induced by the optimal transport. A Gaussian law is determined by its mean vector m and its covariance matrix Σ and is noted N m,Σ . It is well known that the expression of its density f m,Σ is written as:

f m,Σ (x) = 1 (2π) N/2 det(Σ) 1/2 e -1 2 (x-m) t Σ -1 (x-m) .
The space of n-dimensional Gaussian laws, noted G(n), can be identied with the product between the vector space R n and the space of positive denite matrices SPD(n). The space SPD(n) is an open cone included in the vector space of symmetrical matrices Sym(n). Let {e m i } i<n be the canonical basis of R n and {E i,j } be the canonical basis of M (n). Let

E i,i = E i,i , and E i,j = 1 √ 2 (E i,j + E j,i ), ifi < j.
Using this parametrization, E i,j i≤j forms a basis of the space of symmetric matrices.

The space Sym(n

) is identied to R n(n+1) 2 through the coordinates in the basis E i,j i≤j . The union {x l } = {e m i } i<n ∪ E j,k j≤k is a basis of R n × Sym(n). Let consider φ : R n × SPD(n) → G(n) (m, Σ) → N m,Σ Since R n × SPD(n)
is an open subset of the vector space R n × Sym(n) the situation is identical to the one described in Chapter 4.

The diagonalization Σ = RDR t of Σ induces another parametrization of G(n):

ψ : R n × R n + × SO(n) → G(n) (m, λ, R) → N m,(RD λ R t )
where

D λ =    λ 1 . . . λ n   .
In a certain way, (λ, R) are similar to the polar coordinates, with the dierence that R is not necessary of determinant 1 in polar coordinates. The parametrization ψ is not injective, i.e., given

Σ ∈ SPD(n) there exists several (λ, R) ∈ R n + × SO(n) such that Σ = RD λ R t . If this is not a problem in most of theoretical works,
it is when one comes to applications. Indeed it is important to always use the same representation of each object.

Let G(n) be the set of Gaussian laws N m,Σ such that Σ has distinct eigenvalues. Note that the dierence between G(n) and G(n) has a null measure for standard measures. Let dene

E = {λ ∈ R n + | ∀0 ≤ i < j ≤ n, λ i < λ j },
and

H = {D λ | λ ∈ {+1, -1} n , det(D λ ) = 1}.
. The matrix exponential e (.) denes an homeomorphism from a neighborhood of the null matrix in the space of anti-symmetric matrices and a neighborhood of I in D ∈ SO(n).

Let μψ be the Lebesgue measure on R n × R n + × ASym, where ASym is the set of antisymmetric matrices. The computation of .) and α. Thus, knowing the volume change induced by α and Id 2n ×e (.) at (m, λ, 0) is enough to determined the volume change induced by α. We can normalize the Haar measure on rotations such that Id 2n × e (.) does not introduce local at (m, λ, 0). Then the volume change induced by α and α are equal and given by the Jacobian determinant of α.

dµ φ dα * (µ ψ ) (m, D, I) is made through the following diagram: (R n × R n + × ASym(n), μψ ) (R n × R n + × D, µ ψ) (R n × SPD(n), µ φ ) Id 2n × e (.) α α where α(m, λ, A) = (m, e A D λ (e A ) t ). The local volume change induced by α at (m, λ, 0) ∈ R n × R n + × ASym(n) is a product between local volume change of Id 2n × e (
Let A i,j = 1 √ 2 (E i,j -E j,i ). The partial derivatives of α at (m, λ, 0) in basis B = {e m i } ∪ {e λ i } ∪ {A i,j } and B = {e m i } ∪ E i,j ,
can easily be computed and are given by

∂ α ∂e m i (m, λ, 0) = e m i , ∂ α ∂e λ i (m, λ, 0) = E i,i , ∂ α ∂A i,j (m, λ, 0) = (λ j -λ i )E i,j . (6.8)
Since the E i,j form an orthonormal system, we have:

dµ φ dα * (µ ψ) (m, D λ , I) = i<j≤n |λ i -λ j |, (6.9)
and thus using independency of rotation:

dµ φ dα * (µ ψ) (m, D λ , R) = i<j≤n |λ i -λ j |, (6.10)
when eigenvalues are distinct. This result can be found in Example 8.7, page 158, of [START_REF] Wijsman | Invariant measures on groups and their use in statistics[END_REF].

calculate the dierent scalar product between the vectors of the canonical basis. More precisely, ∀i, j, k:

< e m i , E j,k > F isher N m,D = 0, ∀i, j: < e m i , e m j > F isher N m,D = δ i,j λ -1 i ,
where δ i,j = 1 if i = j and δ i,j = 0 if i = j. In addition, ∀i, j:

E i,i , E j,j F isher N m,D = δ i,j λ -2 i , ∀(i = j, k = l): E i,j , E k,l F isher N m,D = 1 2 δ i,k δ j,l (λ i λ j ) -1 , ∀j = k: E i,i , E j,k F isher N m,D = 0.
Thus, one can write

e m i E i,i E i,j G F isher N m,D =              1 λ 1 . . . 1 λ 2 1 λ 2 1 . . . 1 λ 2 n 1 2λ i λ j              and now, from G F isher N m,D , we have dvol F isher dµ φ = det(G F isher N m,Σ ) = 1 2 n(n-1) 2 i 1 λ n+2 i = 1 2 n(n-1) 2 det(Σ) -(n+2)/2 . (6.12)
As announced in section 6.4.1, the measure induced by the Fisher metric is a product between a measure on the mean, a measure on eigenvalues and a measure on rotations, i.e.,

dvol F isher = det(G F isher N m,Σ ) dµ φ dα * (µ ψ) dµ ψ = 1 2 n(n-1) 2 i<j≤n |λ i -λ j | i 1 λ n+2 i dLeb R n dLeb R n + dHaar SO(n) .
(6.13)

The Wasserstein metric

The 2-Wasserstein metric is a product metric between the space of means and the space of covariance matrices:

d W 2 (N m 1 ,Σ 1 , N m 2 ,Σ 2 ) 2 = d m (m 1 , m 2 ) 2 + d Σ (Σ 1 , Σ 2 ) 2 , with d m (m 1 , m 2 ) = |m 1 -m 2 |, d Σ (Σ 1 , Σ 2 ) = tr(Σ 1 ) + tr(Σ 2 ) -2tr( Σ 1/2 1 Σ 2 Σ 1/2
1 ). [START_REF] Takatsu | Wasserstein geometry of Gaussian measures[END_REF] proved that this distance is induced by a Riemannian metric. Let vol W be the measure associated with such Riemannian metric.

The Riemannian volume element

The calculation of the volume measure is similar to the calculation of the volume measure associated with the Fisher metric. Let u and v be vectors of R n Sym(n). Let index m and Σ denote the components associated to R n and Sym(n) respectively. The scalar product at (m, Σ) according to [START_REF] Takatsu | Wasserstein geometry of Gaussian measures[END_REF] is given by

< u, v > W m,Σ = u t m v m + tr(ũ Σ Σṽ Σ ), with u Σ = Σũ Σ + ũΣ Σ, v Σ = Σṽ Σ + ṽΣ Σ. (6.14) Let G W N m,Σ
be the matrix of the scalar product in the canonical basis of R n Sym(n).

The density of the Wasserstein volume measure with respect to the Lebesgue measure of

R n × SPD(n) is given by dvol W dµ φ = det(G W N m,Σ ).
The invariance of the Wasserstein metric under the action of rotations implies the invariance of the measure, i.e.,

det(G

W N m,RDR t ) = det(G W N m,D ).
Thus the volume change only need to be computed for the case diagonal matrices. Let D ∈ SPD(n) be a diagonal matrix of eigenvalues λ i . Equations. (6.14) can now be solved.

Since, equations (6.14) are equivalent to

2λ i E i,i = DE i,i + E i,i D, (λ i + λ j )E i,j = DE i,j + E i,j D,
u Σ = P (ũ Σ ), v Σ = P (ṽ Σ ),
where P is a linear operator of Sym(n) whose matrix form is

P : E i,i E i,j      2λ 1 . . . 2λ n (λ i + λ j )     
The dierent terms of the metric are given by:

∀(i, j ≤ k): e m i , E j,k W N m,D = 0, ∀(i, j): e m i , e m j W N m,D = 1, ∀(i, j): E i,i , E j,j W N m,D = δ i,j 1 4λ i ,
where again δ i,j = 1 if i = j and δ i,j = 0 if i = j, ∀(i, j = k):

E i,i , E j,k W N m,D = 0, ∀i < j, k < l: E i,j , E k,l W N m,D = δ i,k δ j,l 1 2(λ i + λ j )
.

Thus, we obtain

e m i E i,i E i,j G W N m,D =             1 . . . 1 1 4λ 1 . . . 1 4λn 1 2(λ i +λ j )            
As for the Fisher measure, we compute

dvol W dµ φ = det(G W N m,Σ ) = 1 2 1 2 n 2 + 3 2 n i 1 λ i i<j 1 (λ i + λ j ) , (6.15)
and the measure induced by Wasserstein metric is written as

dvol W = det(G W N m,Σ ) dµ φ dα * (µ ψ) dµ ψ = 1 2 1 2 n 2 + 3 2 n i 1 √ λ i i<j |λ i -λ j | (λ i + λ j ) dLeb R n dLeb R n + dHaar SO(n) .
(6.16) 6.4.3.2. The volume change of the exponential map

We are interested here in the computation of the following quantity:

θ Σ 1 : Σ 2 → θ Σ 1 (Σ 2 ) = dvol W dexp * (Leb Σ 1 ) (Σ 2 ),
that is to say the density of the Riemannian measure of G(n) with respect to the Lebesgue measure of the parametrization of G(n) induced by exp Σ 1 . As we have stated, Wasserstein metric is a product between the Euclidean metric on R n and a Riemannian metric on SPD(n). The Euclidean part can be omitted in this section, since it does not aect the volume change.

It has been shown in [START_REF] Takatsu | Wasserstein geometry of Gaussian measures[END_REF] that the L 2 -Wasserstein distance on centered Gaussian measure is induced by a Riemannian metric on SPD(n). Furthermore, the application Π :

Gl n → SPD(n) G → GG t
is a Riemannian submersion when GL n is endowed with the scalar product A, B = tr(AB t ) and SPD(n) is embedded with the Riemannian metric. For A ∈ Gl(n) the kernel of dΠ is called the vertical space, and its orthogonal the horizontal space. For Σ = RDR

t ∈ SPD(n), let Σ 1/2 = RD ( √ λ i ) R t , such that Π(Σ 1/2 ) = Σ.
The horizontal space at Σ 1/2 , see [START_REF] Takatsu | Wasserstein geometry of Gaussian measures[END_REF], is given by:

H Σ 1/2 = {XΣ 1/2 , X ∈ Sym}.
H Σ 1/2 can be identied with T Σ M via dΠ. Let Π A (Z) = Π(A + Z). For simplicity reasons, the exponential map will be expressed in H Σ 1/2 . For Z ∈ H Σ 1/2 classical results on Riemannian submersions give that the image by Π of the geodesic Σ 1/2 + tZ is the geodesic exp Σ (tZ). Hence, one has exp Σ (Z) = Π Σ 1/2 (Z).

dvol W dexp * (Leb H Σ 1/2 ) = dvol W dLeb SPD(n) dLeb SPD(n) dexp * (Leb H Σ 1/2 ) ,
where

dvol W dLeb SPD(n)
was computed in Eq. (6.15). The computation of

dLeb SPD(n) dexp * (Leb H Σ 1/2 ) can
be made through the Jacobian of the exponential application. The computation of the dierential of ΠA gives:

d ΠA (XA, X ∈ Sym(n))(Y A, Y ∈ Sym(n)) = (A + XA)(Y A) t + Y A(A + XA) t = AA t Y + XAA t Y + Y AA t + Y AA t X = AA t Y + Y AA t + XAA t Y + Y AA t X.
The invariance of the Wasserstein metric under the action of rotations enables to restrict the computation of θ Σ 1 to θ D with D a diagonal matrix. So,

d ΠD 1/2 (XD 1/2 )(Y D 1/2 ) = DY + Y D + XDY + (XDY ) t = (I + X)DY + Y ((I + X)D) t .
Y → d ΠD 1/2 (XD 1/2 )(Y D 1/2 ).
Theorem 6.4.2 Let A be a square matrix of order n and

S A : Sym(n) → Sym(n) Y → Y A + A t Y Then det(S A ) = i j (λ i + λ j ) = 2 n det(A) i<j (λ i + λ j ),
where the λ i are the complex eigenvalues of A. Furthermore, det(S A ) depends only on the coecients of the characteristic polynomial of

A χ A (λ) = det(λI -A) = n i=1 (λ -λ i ) = λ n + n k=1 (-1) k σ n-k λ k . 1. If n = 2, det(S A ) = 4(detA)(traceA); 2. If n = 3, det(S A ) = 2 3 detA(σ 1 σ 2 -σ 3 ); 3. If n = 4, det(S A ) = 2 4 detA(σ 1 σ 2 -σ 3 -σ 4 σ 2 1 -σ 2 3 );
4. If n ≥ 4, det(S A ) can be obtained with the help of standard libraries of mathematical programming languages using the decomposition on the elementary symmetric polynomials.

Proof 1. Suppose A is diagonalisable and (X i ) 1≤i≤n is a basis of eigenvectors of A.

Then S A (X i X t j + X j X t i ) = (λ i + λ j )(X i X t j + X j X t i ),
and

(X i X t j + X j X t i ) 1 i j n is a basis of Sym(n). Consequently det(S A ) = i j (λ i + λ j ) = 2 n det(A) i<j (λ i + λ j )
2. det(S A ) is a symmetric polynomial of the λ i . The fundamental theorem of symmetric polynomial states that any symmetric polynomial can by expressed as a polynomial function of the elementary symmetric polynomials. Since the coecients of the characteristic polynomial of A are the elementary symmetric polynomials in λ i , det(S A ) is a polynomial function of the coecient of the characteristic polynomial of A.

3. In the general case, any matrix A can be seen as a limit of diagonalizable matrices A n . By continuity we have det(S A ) = lim det(S An ) and the coecients of the polynomial χ A are the limits of the coecients of the polynomials χ An . Therefore the formula can be extended to any matrix.

The basis B = (E i≤j D 1/2 ) of H D 1/2 is orthogonal for the scalar product A, B = tr(AB t ), with E ii D 1/2 , E ii D 1/2 = λ i and E i<j D 1/2 , E i<j D 1/2 = λ i +λ j 2 . Thus in an orthonormal basis of H D 1/2 , we have

det(B) = λ i i<j λ i + λ j 2 .
Let consider the map f :

H D 1/2 → Sym(n) Y D 1/2 → Y Expressed in orthonormal basis, det(f ) = det(B) -1 . Since d ΠD 1/2 (XD 1/2 ) = S (I+X)D • f , the determinant of d ΠD 1/2 expressed in orthonormal basis is given by det(d ΠD 1/2 ) = det(S (I+X)D )det(f ) = det(S (I+X)D ) det(B) .
Hence the expression of the volume change is given by,

dLeb SPD(n) dexp * (Leb H Σ 1/2 ) = det(d ΠD 1/2 ) = det(S (I+X)D ) det(B) .
(6.17)

Density estimation on multivariate Gaussian laws

Given two matrices Σ 0 , Σ 1 ∈ SPD(n), [START_REF] Takatsu | Wasserstein geometry of Gaussian measures[END_REF] provides the expression of the vector Z ∈ H Σ 1/2 0 :

Z = (Σ 1/2 1 (Σ 1/2 1 Σ 0 Σ 1/2 1 ) -1/2 Σ 1/2 1 -I)Σ 1/2 0 , such that Π Σ 1/2 0 (Z) = Σ 1 . Thus, θ D (Σ) = dvol W dLeb SPD(n) dLeb SPD(n) dexp * (Leb H D 1/2 ) = 1 2 n(n+1) det(S A ) λ Σ i i<j (λ Σ i + λ Σ j ) λ D i i<j (λ D i + λ D j )
,

where (λ D i ) and (λ Σ i ) are the eigenvalues of D and Σ respectively, and where

A = Σ 1/2 (Σ 1/2 DΣ 1/2 ) -1/2 Σ 1/2 D 1/2 .
We recall that A corresponds to the matrix (I + X)D of Eq. (6.17).

Finally, using the invariance of the Wasserstein metric under the action of rotations, we have for Σ 0 = RDR t :

θ Σ 0 (Σ 1 ) = 1 2 n(n+1) det(S A ) λ Σ 1 i i<j (λ Σ 1 i + λ Σ 1 j ) λ Σ 0 i i<j (λ Σ 0 i + λ Σ 0 j ) (6.18)
where (λ Σ 0 i ) and (λ Σ 1 i ) are the eigenvalues of Σ 0 and R t Σ 1 R respectively,

A = (R t Σ 1 R) 1/2 ((R t Σ 1 R) 1/2 D(R t Σ 1 R) 1/2 ) -1/2 (R t Σ 1 R) 1/2 D 1/2 ,
and det(S A ) is given in Theorem 6.4.2.

Density estimation on multivariate Gaussian laws

Let Ω be a space endowed of with a σ-algebra and a probability measure p. Let X be a random variable Ω → G(n). We assume that the support of X is a compact set, noted Supp. The space of Gaussian laws G(n) is equipped successively with the Fisher and the Wasserstein metric. The Riemannian measure is called vol and the measure on G(n) induced by X is noted µ X . We assume that µ X has a density, noted f , with respect to vol, and that the support of X is a compact set noted Supp. Let (x 1 , .., x k ) ∈ M k be a set of draw of X and µ k = 1 k i δ x i denotes the empirical measure of the set of draws. To each point x ∈ G(n) quantities m x and Σ x representing the associated mean and covariance are associated.

The Fisher metric

The expression of distances is unknown, there are no known interesting tilings, and there is no expression of the eigenfunctions of the Laplace operator. Thus there is no way of estimating densities based on a fully Riemannian approach. The only option here is to perform a density estimation in an Euclidean context, and to adapt the result to the Riemannian measure. If X is seen as a random variable in SPD(n) the estimation gives:

f F isher = f Eucl dµ φ dvol = 2 n(n-1) 2 i λ n+2 i f Eucl , (6.19)
where f Eucl is an estimation of the density of X with respect to the measure induced by the scalar product A, B = tr(AB t ). However the estimation does not respect homogeneity and isotropy considerations.

The Wasserstein metric

As with the Fisher metric, there are no known interesting tiling, and there is no expression of the eigenfunctions of the Laplace operator. However it is possible to obtain expressions of kernels. Given a scaling parameter r > 0, the estimator becomes: (6.20) where θ Σx i (Σ x ) was computed in Eq. ( 6.18).

f r k (x) = 1 k i 1 r n 1 θ Σx i (Σ x ) K d W 2 (x, x i ) r ,

The space of multivariate centred Gaussian laws

We consider specically the study of multivariate Gaussian laws, after adding the property m = 0. The set of multivariate centered Gaussian laws is noted G m=0 (n). In G m=0 (n), a Gaussian laws is determined by its covariance matrix Σ only. It elements are noted N Σ . The expression of its density f Σ is:

f m,Σ (x) = 1 (2π) N/2 det(Σ) 1/2 e -1 2 (x) t Σ -1 (x) .
Everything written in section 6.4 remains valid after removal of the mean. The measures µ φ and µ ψ become respectively:

µ φ = Leb SPD(n) , µ ψ = Leb R n + × Haar SO(n)
, and G m=0 (n) is identied with SPD(n), such that N Σ and Σ are used indierently.

The dierent expressions associated to the Wasserstein metric computed in section 6.4.3 do not change for the case of centered Gaussian laws. Furthermore, space G(n) does not gain any property after imposing m = 0 in the Wasserstein metric. Thus this section focuses on the Fisher metric.

The Fisher metric on G m=0 (n) happens to be the metric induced by the Fisher metric on G(n). On the contrary of G(n), G m=0 (n) is symmetric and is then easier to study. Let GL(n) be the set of invertible matrices of size n. For any G in GL(n), the application Σ → GΣG t is an isometry of G m=0 (n). This application corresponds to a coordinate change for covariance matrices. The action of GL(n) on G m=0 (n) induced by this application is transitive. Indeed, for any Σ ∈ SPD(n), let Σ 1/2 ∈ GL(n). The orbit of the identity I is SPD(n): Σ = Σ 1/2 IΣ 1/2 . The stabilizator of I are the orthogonal matrices. Thus G m=0 (n) can be identied to GL(n)/O(n).

There is an explicit expression of the distance between two laws:

d F isher (N Σ 1 , N Σ 2 ) = ||log(Σ -1/2 1 Σ 2 Σ -1/2 1 )||,
where ||.|| is the norm associated with the scalar product tr(AB t ).

The Riemannian volume element

By adapting the calculations done in section 6.4.2, we obtain that the metric G F isher

N D ex- pressed in the canonical basis of SPD(n) is given by E i,i E i,j G F isher N D =       1 λ 2 1 . . . 1 λ 2 n 1 2λ i λ j       Therefore, using the classical machinery, dvol F isher dµ φ = det(G F isher N Σ ) = 1 2 n(n-1) 2 i 1 λ n+1 i = 1 2 n(n-1) 2 det(Σ) -(n+1)/2 . (6.21)
This expression can be found in [START_REF] Terras | Harmonic Analysis on Symmetric Spaces and Applications II[END_REF].

Furthermore, the Riemannian volume element is now

dvol F isher = det(G F isher N Σ ) dµ φ dα * (µ ψ) dµ ψ = 1 2 n(n-1) 2 i<j≤n |λ i -λ j | i 1 λ n+1 i dLeb R n + dHaar SO(n) .
(6.22) 6.6.2. The volume change of the exponential map

We are interested here in the computation of the following quantity:

θ Σ 1 : Σ 2 → θ Σ 1 (Σ 2 ) = dvol F isher dexp * Σ 1 (Leb Sym(n) ) (Σ 2 ),
that is to say the density of the Riemannian measure of G m=0 (n) with respect to the Lebesgue measure of the parametrization induced by exp Σ . Since SPD(n) is homogeneous, the volume change can be computed at the identity matrix I. For the ane invariant metric we have that exp I is the classical exponential application on matrices, see [START_REF] Lenglet | Statistics on the Manifold of Multivariate Normal Distributions: Theory and Application to Diusion Tensor MRI Processing[END_REF]. The density is split into two terms:

dvol F isher dexp * (Leb Sym ) = dvol F isher dLeb SPD(n) dLeb SPD(n) dexp * (Leb Sym(n) )
,

where

dvol F isher dLeb SPD(n)
was already computed in Eq. (6.21). Recall the function dened in section 6.4:

α : R n × SO(n) → Sym(n) (λ, R) → RD λ R t
Let ẽ xp be dened as:

ẽ xp : R n × SO(n) → R n + × SO(n) (λ, R) → e λ i , R
The computation of

dLeb SPD(n) dexp * (Leb Sym(n) ) is made using the following commutation diagram, (R n × Q, µ ψ = Leb(R n ) × Harr SO(n) ) (R n × Q, µ ψ) (Sym(n), µ φ = Leb(Sym)) (SPD(n), µ φ ) ẽ xp α exp α that is to say exp•α = α • ẽ xp.
For the sake of simplicity, µ φ is here extended to Sym(n). Knowing the volume changes induced by α and ẽ xp enables the computation of the volume change induced by exp. The local volume change induced by α was computed in section 6.4 and gives:

dµ φ dα * (µ ψ) = i<j |λ j -λ i |.
We address now the volume change

dexp * (µ ψ ) dµ ψ
induced by ẽ xp, where exp * (µ ψ) denotes the image measure of µ ψ by ẽ xp.

Let F : R n × Q → R, F (λ, R) = f (λ)g(R) such that R n ×Q F (λ, R)dexp * (µ ψ) = R n ×Q F ( ẽ xp(λ, R))dµ ψ = R n ×Q F (e λ , R)dµ ψ = R n f (e λ )dλ Q g(R)dHaar SO(n) = R n f (Λ) n i=1 1 Λ i dΛ Q g(R)dHaar SO(n) = R n ×Q F (Λ, R) n i=1 1 Λ i dµ ψ.
Since the set of separable functions F ∈ L 1 is dense in L 1 , we have that:

dexp * (µ ψ) dµ ψ (λ, R) = 1 λ i and dµ ψ dexp * (µ ψ) (λ, R) = λ i . Using dµ φ dexp * (µ φ ) (RD e λ R t ) dµ φ dµ ψ (λ, R) = dµ φ dµ ψ (e λ , R) dµ ψ dexp * (µ ψ) (e λ , R),
we have that

dµ φ dexp * (µ φ ) (RD e λ R t ) i<j |λ i -λ j | = i<j |e λ i -e λ j |e i λ i ,
and when eigenvalues are distinct,

dµ φ dexp * (µ φ ) (RD e λ R t ) = dLeb SPD(n) dexp * (Leb Sym(n) ) (RD e λ R t ) = i<j |e λ i -e λ j | |λ i -λ j | e i λ i . (6.23)
The continuity of the volume change enables to extend the expression everywhere. Let log(Σ) = RD (log(λ i )) R t , where again Σ = RDR t ∈ SPD(n). The application log is the inverse the exponential application on SPD(n). Combining Eq. ( 6.23) and Eq. ( 6.21), we obtain for Σ ∈ SPD(n):

θ I (Σ) = 1 2 n(n-1)/2 1 e (n+1)λ i |e λ j -e λ i | |λ j -λ i | e λ i = 1 2 n(n-1)/2 |e λ j -e λ i | |λ j -λ i | e (1-n+1 2 )λ i = 1 2 n(n-1)/2 1 |λ i -λ j | 2e λ i +λ j 2 |e λ i -λ j 2 -e λ j -λ i 2 | 2 e (1-n+1 2 )λ i = 2 n(n-1)/2 2 n(n-1)/2 1 |λ i -λ j | sinh |λ i -λ j | 2 e λ i +λ j 2 e (1-n+1 2 )λ i = 1 |λ i -λ j | sinh |λ i -λ j | 2 e (n-1) 2 λ i e (1-n+1 2 )λ i = sinh |λ i -λ j | 2 |λ i -λ j | ,
where the λ i are the eigenvalues of log(Σ). This expression is not new since the factor in sinh(.) is known to be the Riemannian volume element in polar coordinates. Let Σ 1/2 = RD ( √ λ i ) R t . Using the invariance of the metric we obtain nally, , (6.24) where the λ i are the eigenvalues of log(Σ

θ Σ 1 (Σ 2 ) = 2 n(n-1)/2 sinh |λ i -λ j | 2 |λ i -λ j |
1/2 1 Σ 2 Σ 1/2
1 ).

FourierHelgason transform

The G m=0 (n) is symmetric and there is a FourierHelgason transform dened on it. This paragraph is mainly extracted from Peter T. [START_REF] Kim | Deconvolution Density Estimation on the Space of Positive Denite Symmetric Matrices[END_REF]. Let

p s : SPD(n) → C Σ → i |Σ| s i i (6.25)
where |Σ| i is the determinant of Σ after removing row and column i.

6.7. Density estimation on multivariate centred Gaussian laws 6.7.1. Histograms

Let (A j ) be a partition of G m=0 (n). Let α j be the number of draws that lay in A j , α j = Card{x i ∈ A j }. The histogram of the draws (x i ) is the function

fk = 1 k j 1 vol(A j ) α j 1 A j ,
where 1 E is the indicator function of E. Thus fk is an estimator of the density f . In order to respect the criterion mentioned in section 6.2, bins should be isometric and verify some isotropy properties. In the Euclidean context, square bins present an interesting compromise between isotropy and computational simplicity. The main drawback of this method is the diculty to obtain a regular partition of the space. A lead is to see G m=0 (n)

as the imaginary part of the Siegel upper half space. See [START_REF] Zelikin | Control Theory and Optimization I Homogeneous Spaces and the Riccati Equation in the Calculus of Variations[END_REF] for an introduction to the Siegel upper half space. Regular polygonal tiling on the hyperbolic space, which can be extended to the Siegel space, are widely studied and are good candidates for the histogram bins. However the absence of homothetic transforms makes the adaptation of the estimation to the scaling factor dicult.

Kernel density estimation

The explicit expression of the change of volume θ p (q) induced by the exponential map enables us to use the kernel density estimation. Given a scaling parameter r > 0, the estimation becomes: (6.30) where θ Σx i (Σ x ) is given in Eq. (6.24). Thus, we have the following expression for the Pelletier kernel density estimator (6.31) where the λ . are the eigenvalues of log(Σ

f r k = 1 k i 1 r n 1 θ Σx i (Σ x ) K d F isher (x, x i ) r ,
f r k = 1 k i 2 n(n-1)/2 r n p<q sinh( |λp-λq| 2 ) |λ p -λ q | K ||log(Σ -1/2 x i Σ x Σ -1/2 x i )|| r ,
-1/2 x i Σ x Σ -1/2 x i
).

Orthogonal series density estimation

As discussed in previous section, the symmetry of the space gives rise to a FourierHelgason transform and its inverse. The empirical transform is given by:

H emp (f )(s, R) = i f (Σ x i )e s,R (Σ x i ),
where e s,R is the plane wave dened in Eq. (6.26).

For s ∈ C n (ρ), see Eq. (6.28), let

r j = s j + s j+1 + • • • + s n + 1 4 (n -2j + 1) and λ s = -(r 2 1 + r 2 2 + • • • + r 2 n - 1 48 n(n 2 -1)), then ∆e s,R (Σ) = λ s e s,R (

Σ).

Let C n (ρ, T ) = {s ∈ C n (ρ), λ s < T }. For a cut-o parameter T > 0, the associated density estimator is then given by:

f T k (Σ) = s∈C n (ρ,T ) R∈Q H emp (f )(s, R)e s,R (Σ)dHaar SO(n) dµ s . (6.32) Let 1 F T = C n (ρ, T ) × SO(n).
It can be veried that H -1 (1 F T ) is invariant under the action of SO(n). Thus, one has

f T k = µ k * H -1 (1 F T ).
In other words, the estimation do not privilege specic directions or locations. The main drawback of this method is its computational cost: one has to evaluate integral (6.32). The integral present in the transform is due to the non-compactness of the underlying space SPD(n). In R n , the assumption that the random variable lays in a parallelepiped P , that is to say Supp ⊂ P , enables the use of Fourier series instead of the Fourier transform. Then the estimation involves only a nite sum. Unfortunately, there are no compact domains D ⊂ SPD(n) where the eigenfunctions of the Laplacian are known.

The space of Gaussian laws with xed rotation

Given a rotation matrix R, we study here the specicity of the space of Gaussian laws of the form N (m,RDR t ) where D is diagonal matrix, under the Wasserstein metric. We choose here to represent the Gaussian law N (m,RDR t ) by its mean m ∈ R n and a set of standard deviation σ ∈ R n , with σ 2 i an eigenvalue of D. In the (m, σ) parametrization, the expression of the metric found in (6.4.3.1) takes the following form:

e m i e σ i G W asserstein Nm,σ =           1 . . . 1 1 . . . 1          
Thus in the (m, σ) parametrization, the metric is Euclidean. The density estimation can then be achieved using standard techniques.

The space of Gaussian laws with xed mean and rotation

Given a mean m and a rotation R, we study here the specicity of the space of Gaussian laws of the form G m,R (n) = {N (m,RDR t ) , with D diagonal}, under the Fisher metric. The metric computed in (6.4.2) expressed in the basis {e λ i } gives:

e λ i G F isher N m,RD λ R t =    1 λ 2 1 . 1 λ 2 n   
After the change of coordinates γ i = log(λ i ), the metric becomes:

6.10. The space of Gaussian laws with xed covariance

e γ i G F isher N m,RD λ R t =   1 . 1  
Thus in the γ = log(λ) parametrization, the metric is Euclidean. The density estimation can then be achieved using standard techniques.

6.10. The space of Gaussian laws with xed covariance Given a covariance matrix Σ, we study now the specicity of the space of Gaussian laws of the form N (m,Σ) , under the Fisher metric. Note that the mean under the Wasserstein metric follows the canonical Euclidean scalar product. From Eq. (6.11) in G(n):

< u, v > F isher m,Σ = u t m Σ -1 v m + 1 2 tr(Σ -1 u Σ Σ -1 v Σ ),
and we have that, at a xed Σ, it is given by

< u m , v m > F isher m = u t m Σ -1 v m .
The metric on the mean is Euclidean, of scalar product Σ -1 . The density estimation can then be achieved using standard techniques.

The space of univariate Gaussian laws

This section addresses the case of G(n = 1). Space G(n = 1) is identied with the half plane R × R + . That means that each normal law N is described by its mean m ∈ R and its (positive) standard deviation σ ∈ R + . Note that, as in section 6.4.3, there is a slight change of convention with respect to the previous sections since the equivalent of the covariance matrix Σ would be the square of the standard deviation σ 2 . The study is focused on the Fisher metric since the case of the Wasserstein metric has been analyzed in section 6.8.

Fisher metric and the Poincaré upper half plane

The Poincaré upper half plane of curvature R is a half plane {x ∈ R, y ∈ R + } endowed with the following Riemannian metric:

G (x,y) = a 2 1/y 2 0 0 1/y 2 , (6.33) with R = -1 a 2 .
The Poincaré upper half plane is a model of hyperbolic geometry, see Can- non et al. (1997). The distance between (x 1 , y 1 ) and (x 2 , y 2 ) is given by: d((x 1 , y 1 ), (x 2 , y 2 )) = a cosh -1 1 + (x 2 -x 1 ) 2 + (y 2 -y 1 ) 2 2y 1 y 2 .

(6.34)

The Poincaré upper half plane is isotropic: each location and directions are equivalent.

After having considered the change of convention between Σ and σ, the matrix of the Fisher metric on G(n = 1) can be deduced from Eq. ( 6 where ds is the Jacobian matrix of s, the expression of the metric is obtained in the new parametrization:

G F isher (x,y) = 2 1/y 2 0 0 1/y 2 (6.37)
Thus, under the Fisher metric, the space G(n = 1) is isometric to a Poincaré upper half plane of curvature -1 2 . Let H denotes just the Poincaré upper half plane of curvature -1 2 .

The Riemannian volume element

In the (x, y) parametrization, the Riemannian volume measure vol has the following expression:

dvol dLeb (x,y) = det(G F isher (x,y) ) = 2 y 2 .
(6.38) 6.11.3. The volume change of the exponential map Given a reference point p, any point in polar coordinates (r, α) of the hyperbolic space is dened as the point q (r,α) at distance r of p on the geodesic with initial direction α ∈ S 1 .

Since the hyperbolic space is isotropic, the expression of the length element in polar coordinates depends only on r, see [START_REF] Anker | The Heat Kernel on Noncompact Symmetric Spaces[END_REF]; [START_REF] Grigoryan | Heat Kernel and Analysis on Manifolds[END_REF]. Expressed in polar coordinates the hyperbolic metric of curvature -1 is:

g = dr 2 + sinh(r) 2 g S 1 .
For a hyperbolic space of curvature -1 2 , the expression becomes: g = 2(dr 2 + sinh(r) 2 g S 1 ).

In the basis e r = dp (r,α)

dr , e α = dp (r,α) rdr , the metric takes the following form:

G (r,α) = 2 1 0 0 sinh(r) 2 1 r 2 (6.39)
Thus, if we consider vol as a measure on the parametrization (r, α), one has 6.40) where r = d(p, q). The polar coordinates are a polar expression of the exponential map at p. Indeed the exponential map at p can be dened as exp p (r(cos(α), sin(α))) = q (r,α) .

dvol dLeb (r,α) (r, α) = det(G (r,α) ) = 2 sinh(r) r , (
6.12. Density estimation on univariate Gaussian laws under the Fisher metric Therefore, considering now vol as a measure on the Riemannian space, the volume change of the exponential is just: dvol dexp * (Leb p ) (q (r,α) ) = 2 sinh(r) r .

(6.41)

We obtain then θ p (q) = 2 sinh(d(p, q)) d(p, q) , (6.42)

where θ is the function used in Eq. ( 6.3).

FourierHelgason transform

We have said that G(n = 1) is symmetric, therefore there exists a FourierHelgason transform. This paragraph is mainly extracted from [START_REF] Huckemann | Möbius deconvolution on the hyperbolic plane with application to impedance density estimation[END_REF]. Space H is identied with the complex upper half plane: (1, 0) = 1 and (0, 1) = i. For (s, k) ∈ C × SO(2), the Helgason functions (plane waves) have in the case the following expression:

e s,ku : H → C z → Im(k u .z) s = e Re(s) log(Im(ku.z)) e iIm(s) log(Im(ku.z)) (6.43) where Re(.), Im(.) denotes the real and imaginary part, k u = cos(u) sin(u) -sin(u) cos(u) , and k u .z = z cos(u) + sin(u) cos(u) -z sin(u) .

We have that ∆e s,ku = s(s -1)e s,ku , see [START_REF] Huckemann | Möbius deconvolution on the hyperbolic plane with application to impedance density estimation[END_REF]. Let f ∈ C ∞ c (H), where C ∞ c (H) denotes the set of innitely dierentiable with compact support. For (s, k) ∈ C × SO(2), the FourierHelgason transform of f is given by:

FH(f )(s, k u ) = H f (z)e s,ku (z)dvol. (6.44)
The inverse transform is given by:

f (z) = t∈R u=2π u=0 FH(f )( 1 2 + it, k u )e 1 2 +it,ku (z) 1 8π 2 t tanh(πt)dudt.
(6.45) 6.12. Density estimation on univariate Gaussian laws under the Fisher metric Let X be a random variable in Ω → G(n = 1). Notations of this section are similar to notation of section 6.2, except for the set of draw which now is (z 1 , ..., z k ), with z i = (m z i , σ z i ). For the reasons explained in the previous section we are here only interested in the Fisher metric.

other hand, as mentioned in section 6.2.5 the curse of dimensionality pushes us to reduce the dimension of the studied objects.

We start the discussion with an example. Let X be a random variable valued in R 2 of density f . Consider the distribution of the radius of the polar coordinates (r, θ). Firstly, the Lebesgue measure of R 2 can be written as a product dLeb = rdθdr. Thus the Lebesgue measure of R 2 induces a measure on radiuses rdr. The density computed with respect to rdr can be interpreted as an average of f over a slice of constant radius. Secondly, the metric on R 2 induces a natural metric on the space of radius seen as the quotient space R ∼ R 2 /θ. The natural quotient metric is the Euclidean metric on R, the associated measure being the Lebesgue measure dr. Thus this example shows that there is not a unique way of addressing the density estimation of partial quantities.

Remind the context of section 6.2. Let Ω be a space, endowed of with a σ-algebra and a probability measure p. Space G(n) is equipped with a Riemannian metric G, the associated Riemannian measure is called vol. Let X be a random variable, X, Ω → G(n). The measure on G(n) induced by X is noted µ X . We assume that µ X has a density, noted f , with respect to vol. G(n) can be identied to the product R n × E × D, see section 6.4. Recall that vol(G(n) \ G(n)) = 0. Under this identication, the measure vol is described by a product vol = µ mean × µ λ × µ R , see section 6.4.1. Let X mean , X λ and X R be the random variables naturally induced by the identication between G(n) and R n × E × D. Note that X R is not always properly dened. This problem can be neglected given that vol(G(n) \ G(n)) = 0 and that the law of X has a density with respect to vol. Let µ Xmean , µ X λ , µ X R be the associated measures.

Average over slices

Measure µ Xmean (A) is the probability that m x ∈ A when x follows µ X . Thus, µ Xmean (A) = µ X (A × E × D).

Let f mean be the density of µ Xmean with respect to µ mean . Since

µ X (A × E × D) = A E D f dµ R dµ λ dµ mean , we have, f mean (m) = E D f (m, λ, R)dµ R dµ λ . Let α = (µ λ (E)µ R (D)) -1 .
Quantity αf mean can be interpreted as an average value of the original density f over a slice E × D. The same hold for µ X λ and µ X R . Let µ Eucl • be the measure induced by the canonical Euclidean metric of R n × Sym(n). Up to a scaling factor, we recall here the expression of the dierent measures:

µ Eucl mean = Leb, dµ Eucl λ = i<j≤n |λ i -λ j |dLeb(λ), µ Eucl R = Haar SO(n) , µ F isher mean = Leb, µ F isher R = Haar SO(n) , µ W asserstein mean = Leb, dµ W asserstein λ = i 1 √ λ i i<j |λ i -λ j | (λ i + λ j ) dLeb, µ W asserstein R = Haar SO(n) ,
where Haar SO(n) is restricted to D. The expression of µ F isher λ varies between G(n) and G m=0 (n). When G(n) is the underlying space, we have: The question of the underlying metric associated to these measures has not yet been claried. In the absence of underlying metric, the notions of homogeneity and isotropy disappear. Thus the estimation can be achieved using the approach described section 6.2.1.

dµ F isher λ = i<j≤n |λ i -λ j | i 1 λ n+2 i dLeb,
Despite the absence of clear argument, the Lebesgue measure and the Haar measure push to the use of an underlying translation invariant metric.

Rotation

Given a Gaussian law N m,λ,R , the Euclidean metric of R n × Sym(n), the Fisher and the Wasserstein metric induce a left-invariant metric on Q:

d N m,λ,R (R 1 ∈ Q, R 2 ∈ Q) = d(R 1 . Lef t N m,λ,R , R 2 . Lef t N m,λ,R ).
(6.53)

However according to the expressions of the metrics in the basis (E i≤j ) and the dierential of α, see Eq. 6.8, it can be veried that the left-invariante metric is dependent of the choice of N m,λ,R . Thus neither the Euclidean metric of R n × Sym(n), nor the Fisher nor the Wasserstein metric induce a canonical quotient metric on SO f (z) = d F isher (z, (0, 1)) 2 -11 d F isher (z,(0,1))<1 .

(n) ∼ G(n)/(R n × R n + ).
From the draws, rstly, the density is estimated using an Euclidean kernel method, followed by an adaptation to the Riemannian measure, as described in section 6.2.1. Secondly, the density is estimated using a Riemannian kernel, see section 6.2.3. The base kernel is a

quartic kernel K(||x||) = 3 π (1 -||x|| 2 ) 2 1 ||x||<1 .

Density estimation in structure tensor images

The structure tensor image from a grayscale image is a eld of symmetric positive denite matrices. Despite the fact that these matrices are not directly covariance matrices of Gaussian laws, they can be studied under the Fisher and the Wasserstein metric of G m=0 (n), typically n = 2 for 2D images and n = 3 for 3D ones. Fig. 6.7(b)-(g) shows the densities of the two eigenvalues λ 1 ≤ λ 2 of the structure tensor eld computed from Fig. 6.7(a).

For each metric, the densities of eigenvalues are computed following considerations of section 6.13. For the induced measures µ Eucl λ , µ F isher λ and µ W asserstein λ , the densities are rstly estimated using an Euclidean kernel density estimation, followed by an adaptation to the reference measures, given in section 6.13.1. According to section 6.13.2.2, each quotient metric is Euclidean in the right parametrization. In the adapted parametrizations, the densities are obtained using an Euclidean kernel density estimation. The limiting factor for the use of density estimation is the poor spatial resolution. Indeed, the spatial resolution is not yet sucient to enable the estimation of densities on local spatial neighborhoods.

Radar density estimation on the Poincaré disk

In this section the random variable is not an image I but a set of coecients derived from a radar signal. In particular, the draws (x i ) are not explicitly Gaussian laws, but falls into the Poincaré disk, an analog of the Poincaré upper half plane.

Let us discuss briey how radar data are related to hyperbolic space via reection coefcients, for more details see [START_REF] Arnaudon | Riemannian Medians and Means With Applications to Radar Signal Processing[END_REF]. Each radar cell is a complex vector z = (z 0 , • • • , z n-1 ) considered as a realization of a centered stationary Gaussian process

Z = (Z 0 , • • • , Z n-1 ) of covariance matrix R n = E[ZZ * ]. The matrix R n has a Toeplitz structure. For 1 ≤ k ≤ l ≤ n -1, the k-th order autoregressive estimate of Z l is given by Ẑl = -k j=1 a (k) j Z l-j , where the autoregressive coecients a (k) 1 • • • a (k) k are chosen such
that the mean squared error E(|Z l -Ẑl | 2 ) is minimized. In practice, reection coecients are estimated by regularized Burg algorithm BARBARESCO and [START_REF] Barbaresco | Super Resolution Spectrum Analysis Regularization: Burg, Capon and Ago-antagonistic Algorithms[END_REF]. The last autoregressive coecient a

k is called the k-th reection coecient, denoted by µ k and which has the property |µ k | < 1. The coecient for k = 0 corresponds to the power, denoted P 0 ∈ R * + . The reection coecients induce a (dieomorphic) map ϕ between the Toeplitz Hermitian positive denite (HPD) matrices of order n, T n , and reection coecients:

ϕ : T n → R * + × D n-1 , R n → (P 0 , µ 1 , • • • , µ n-1 )
where D = {ζ ∈ C : |ζ| < 1} is the open unit disk of the complex plane. Dieomorphism ϕ is very closely related to theorems of [START_REF] Trench | An Algorithm for the Inversion of Finite Toeplitz Matrices[END_REF].

The Riemannian geometry of the space of reection coecients has been explored in [START_REF] Barbaresco | Information Geometry of Covariance Matrix: Cartan-Siegel Homogeneous Bounded Domains, Mostow/Berger Fibration and Fréchet Median[END_REF] through the Hessian of Kähler potential. The metric appropriate metric on

T n is ds 2 = n dP 2 0 P 2 0 + n-1 k=1 (n -k) |dµ k | 2 (1 -|µ k | 2 ) 2 .
(6.54)

According to the metric (6.54) the space T n can be seen as a product of the Riemannian manifold R * + , ds 2 0 , with ds 2 0 = n(dP 2 0 /P 2 0 ) (logarithmic metric multiplied by n), and (n -1) copies of D, ds 2 k 1≤k≤n-1 , with ds 2 k = (n -k)ds 2 D . R * + × D n-1 , ds 2 is a Cartan Hadamard manifold whose sectional curvatures are bounded, i.e., -4 ≤ K ≤ 0. This metric is related to information geometry and divergence functions, see discussion in [START_REF] Barbaresco | Koszul Information Geometry and Souriau Geometric Temperature/Capacity of Lie Group Thermodynamics[END_REF]. From the product metric, closed forms of the Riemannian distance, arc-length parameterized geodesic, etc. can be obtained, see [START_REF] Barbaresco | Information Geometry of Covariance Matrix: Cartan-Siegel Homogeneous Bounded Domains, Mostow/Berger Fibration and Fréchet Median[END_REF]; Arnaudon et al. and Meer (2009a). However in order to gain generality authors chose to work with pseudo kernels in the sense that they do not take into account the volume change between the exponential map and the manifold. Indeed, this term is generally unknown. The study can thus be resumed for spaces Gaussian laws, except multivariate Gaussians under the Fisher metric. Especially since this thesis focuses on mathematical morphology and on density estimation, it is important to point out the similarity between the mean shift algorithm and the watershed transform of the density. Indeed the mean shift can be interpreted as a steepest ascent on the density.

Probability densities are also important tools in Bayesian classication. The classication decision consists in maximizing the probability of a class given an observation. This requires a priori information on classes and probability densities knowing the class. Our future research on this topic will study the inuence of the choice of the metric on classication results. This type of classication has been addressed on SPD(n) in a recent work [START_REF] Said | Riemannian Gaussian Distributions on the Space of Symmetric Positive Denite Matrices[END_REF], where they model conditional densities by (parametric) Gaussian distributions on SPD(n). We remind that our approach here allows us a non-parametric density estimation on SPD(n).

Part III. for mathematical morphology: given a total order, there are always images where morphological operators will introduce irregularities and aliasing on processed images. In this thesis, the notion of quality of orders and operators is restricted to topological considerations. Given an image it is possible to asses the quality of a total order regarding its 122 7.1. Conclusions impact on morphological operators. This quality measure is represented by a cost function introduced in Chapter 2. Thus given an image it is natural to search orders optimizing this quality measure. The main drawback of the cost function is its computational complexity. A recursive minimization procedure aiming at minimizing the cost has been proposed.

Results in various contexts have been presented. The orders obtained by our minimization procedure produce interesting morphological erosions and dilations. However, orders have to be computed for each images and the minimization of the cost is a time demanding operation. Thus, the minimization framework does not yet an undeniable interest in front of orders such as the bit-mixing. Indeed the bit-mixing order presents interesting topological properties while being easy to evaluate and independent from the image.

The idea of n-ary morphology introduced in Chapter 3 has emerged as a circumvention of the theoretical limits of total orders: a structure based on only two operators can not properly take into account the geometry of a multidimensional space. On the real line each element can be approached by only two sides. In multidimensional spaces, the topological incompatibility exhibited by Lemma 2.3.1 arises from the increase of directions by which a point can be approached. This corresponds to the connectivity assumption on the complementary of balls. Hence the idea of introducing multiple pairs of erosions and dilations appears naturally. Each pair would then be associated with a pole of the space, each pole playing an equivalent role in the studied image. This underlying idea led to a formalization of n-ary morphology. Recall that n-ary images refer here to images valued in a nite set where no mathematical structure is assumed. In comparison with other works addressing the same type of images, the two main specicities of the present paradigm are the introduction of multiple operators and the absence of an indeterminate label or category. Several elementary properties of the binary framework are recovered. An extension to a continuous case has been proposed. Unfortunately theoretical properties are not recovered as naturally as in the n-ary framework. However, operators produce promising results in practice.

Part II

The formalization of image histograms in terms of density estimation suggests to address the estimation using the coarea formula. This approach has already been proposed recently in [START_REF] Duy | Integrating Isosurface Statistics and Histograms[END_REF], originally aiming at visualization applications.

Perceptual color spaces like the Lab space are built such that their natural Euclidean metric is closer to the perceptual metric than the Euclidean metric of other spaces such as the RGB space. Since density estimation mainly involves local computations, the Riemannian assumption on the perceptual metric on colors suggest to locally adapt computations using the local scalar product. Indeed, as shown in Chapter 5 this enables to be more The minimization procedure of the measure of irregularities of total order is not yet fully convincing. However, we are condent in the fact that the optimization can be signicantly improved. The use of techniques such as evolutionary algorithms might for instance enable to signicantly enhance the number of clusters at each step of the recursion, see [START_REF] Flórez-Revuelta | Topology-Preserving Ordering of the RGB Space with an Evolutionary Algorithm[END_REF]. This is currently an ongoing work.

After having examined the transposition of properties of binary morphology to n-ary morphology, our future works will focus on applications. The main considered areas are material science images and pixelwise classied images. The proposed extension of n-ary morphology to a continuous case is still a preliminary work. Before addressing applications, future researches will continue investigating new continuous operators, in particular openings by reconstruction and levelings.

Part II

Ideas explored in Chapter 5 on images histograms provides several avenues for research:

While more intuitive for functions f : R p → R q≤p the coarea formula can also be formulated of functions f : R p → R q>p . The existing work [START_REF] Duy | Integrating Isosurface Statistics and Histograms[END_REF] only addresses the case of real valued functions. A work can thus be carried on so as to extend their proposition to color and multispectral images.

The preliminary work on perceptual color histograms raises rstly the problem of Riemannian metric interpolation. The main diculty is that despite the primary intuition, there is no intrinsic criteria that enables to evaluate the quality of an interpolation. The quality evaluation of an interpolation necessarily depends on the parametrization. A rst study must then be carried out in order to understand how to determine in which parametrization the interpolation should be performed. We did not nd -at least, not yet -a clear answer to this question in the literature.

Secondly, it would be interesting to study the convergence of the proposed density to the true density of the distribution of colors with respect to the interpolated Riemannian metric.

Chapter 6 describes and gives strategies to estimate densities on spaces of Gaussian laws.

Our future researches will focus on the potential uses of the estimated densities.

The mean shift algorithm proposed by [START_REF] Fukunaga | The estimation of the gradient of a density function, with applications in pattern recognition[END_REF] can be adapted to spaces of Gaussian laws using the expressions of kernels computed in Chapter 6.

Given a point cloud, the mean-shift algorithm is a standard mode seeking technique.

This technique is mainly used for point cloud segmentation tracking and smoothing [START_REF] Subbarao | Nonlinear Mean Shift over Riemannian Manifolds[END_REF]. When used for segmentation purposes, the mean- Abstract: Digital image processing has followed the evolution of electronic and computer science. It is now current to deal with images valued not in {0, 1} or in gray-scale, but in manifolds or probability distributions. This is for instance the case for color images or in diffusion tensor imaging (DTI). Each kind of images has its own algebraic, topological and geometric properties. Thus, existing image processing techniques have to be adapted when applied to new imaging modalities. When dealing with new kind of value spaces, former operators can rarely be used as they are. Even if the underlying notion has still a meaning, a work must be carried out in order to express it in the new context.

The thesis is composed of two independent parts. The first one, "Mathematical morphology on non-standard images", concerns the extension of mathematical morphology to specific cases where the value space of the image does not have a canonical order structure. Chapter 2 formalizes and demonstrates the irregularity issue of total orders in metric spaces. The main results states that for any total order in a multidimensional vector space, there are images for which the morphological dilations and erosions are irregular and inconsistent. Chapter 3 is an attempt to generalize morphology to images valued in a set of unordered labels.

The second part "Probability density estimation on Riemannian spaces" concerns the adaptation of standard density estimation techniques to specific Riemannian manifolds. Chapter 5 is a work on color image histograms under perceptual metrics. The main idea of this chapter consists in computing histograms using local Euclidean approximations of the perceptual metric, and not a global Euclidean approximation as in standard perceptual color spaces. Chapter 6 addresses the problem of non parametric density estimation when data lay in spaces of Gaussian laws. Different techniques are studied, an expression of kernels is provided for the Wasserstein metric.

  Figure 2.2.: (a) Original image, (b) grey-scale representation of the lexicographic order, (c) lexicographic dilation.

Figure 2

 2 Figure 2.3.: Quantitative evaluation of the topological distortion (gure extracted from Chanussot and Lambert (1998)). Five white dots are included in the gray square, the weight assigned to the black point is 5/8. These weights are averaged over all the points of the region of interest, as the RGB cube for color images.

  a and c are two versions of a rst color C 1 . Let b be a third point far from a and c. b represents a second color C 2 , distinct from C 1 . If C I (a, b) = 0 or C I (c, b) = 0, no irregularity is created by the triplet a, b, c regardless of the chosen order. Indeed in the image the color C 2 is in presence of at most one version of the color C 1 . However if C I (a, b) and C I (c, b) are both signicant, the color C 2 is neighboring of the two versions of C 1 . It is then important that a < c < b or b < a < c. Thus the quality of a total order can be measured by evaluating to what extend the following property is respected:

Figure 2

 2 Figure 2.4.: Construction of < 2 (see description of the algorithm).

Fig. 2 .

 2 Fig. 2.5 presents the dierent steps on a simple example.

  Figure 2.5.: This example shows the successive steps of the algorithm launched with threshold = 2. (a) Initial point cloud; (b) depth = 1 ; (c) depth = 2; (d) depth = threshold : points in each cluster are ordered according to their distance to neigbhour clusters.

Figure 2

 2 Figure 2.6.: The function f (x, y) used in the experiment is of the form x • G(y). The parameter r 0 of the gate G is related to the parameter r of the Lemma. Triplet of values a < b < c such that d(a, c) ≥ r 0 has no inuence in the cost function.

  Figure 2.7.: A set of 20 points is ordered according to the lexicographic order in (a), the bit-mixing order in (b), and the adapted order in (c).

  Figure 2.9.: Projection of the total order on the image support: (a) original RGB image I(p); (b) lexicographic order (R → G → B); (c) bit-mixing order; (d) our image adapted total order.

Fig. 2 .

 2 Fig.2.9 represents the RGB color image together with projections of the total orders on the image support for the three studied total orders. As we can see, the blue spots are totally invisible to the lexicographic ordering. Note that the lexicographic ordering starts with red, and ends with blue. The bit-mixing paradigm and the image adapted total order

  Figure 2.10.: Morphological processing of image from Fig. 2.9: openings (top row) γ B (I) and closings (middle row) ϕ B (I) using lexicographic order in (a), bit-mixing order in (b) and our image adapted total order in (c). Structuring element is B is a square of 7×7 pixels. Zoom on aliasing introduced by the lexicographic order in (d) and (e).

  Figure 2.11.: Projection of the total order on the image support: (a) original RGB image I(p); (b) lexicographic order; (c) bit-mixing order; (d) our image adapted total order.

Figure 2 .

 2 Figure 2.15.: Example of DTI image.

  2.16.

Figure 2 .

 2 Figure 2.16.: Standard parametrization of a symmetric matrix.

  the Log-Euclidean metric introduced in Arsigny et al. (2006), i.e., d(A, B) = log(A)log(B) .

Fig. 2 .

 2 Fig.2.17 and Fig.2.18 show respectively the results of openings and closings for the dierent total orders. As in standard morphological processing, opening and closing removes small objects respectively from the image foreground and background.

  Figure 2.20.: Median ltering of a DTI image (reference image in Fig. 2.15 and marker images are the openings in Fig. 2.17): using LEX1 in (a), LEX2 in (b), BM IX1 in (c), BM IX2 in (d), image adapted total order using Frobenius norm in (e), and using Log-Euclidean norm in (f ).

  Figure 2.21.: Time lapse sequence of images from the retina: (a) mean image, (b) standard deviation image, (c) corresponding point cloud in the uni-dimensional Gaussian laws space.

  h (I) denotes the upper level set at height h of function I; i.e., X + h (I) = {p ∈ Ω : I(p) ≥ h} .

  Figure 3.1.: Grey-level morphological processing: (a) original image I, (b) closing ϕ B (I), (c) opening γ B (I).

  Figure 3.2.: (a) A three independent grey-level image and (b) the same image where the grey values has been replaced by colors. (c) Example of multivariate (color) image.

Fig. 3

 3 Fig. 3.2(c). As the red usually has a real signication in terms of a particular class of objects, it is very natural to try to remove the red spot, which is not possible using generic classical morphology. The more the image has a complex semantic structure, such as a color image, the more it is dicult to nd a lattice structure which makes every interesting object an extremum.

  Figure 3.4.: (a) Original image, (b) image after applying the blue erosion: one needs to decide how to ll the gap.

  us set a few notations used in the following. If φ is an operator, let φ k be φ • .... • φ the iteration of φ, k times. Let φ |A be the restriction of φ to the subset A. Let us set

  Figure 3.5.: Opening and closing on a 3-ary image: (a) original image I, (b) opening of red color γ red (I; B), (c) closing of black color ϕ black (I; B).

Fig. 3

 3 Fig.3.6 illustrates the dierence between classical geodesic reconstruction and the proposed n-ary reconstruction. For the classical reconstruction, the 3-label image is simply viewed as a grey-scale image.

  Figure 3.7.: Computing the gradient of a color image: (a) original n-ary image, (b) projection of a total color order on a grey-scale image, (c) classical morphological gradient based on order represented by (b), (d) L 1 norm of the n-ary morphological gradient grad(I)(x) 1 .

  -based erosion: Erosion takes the value of the major label dierent from i in the structuring element B: ε i (I; B)(x) = min(arg max j =i (Card {p ∈ B x |I(p) = j})).

  -based erosion: Erosion replaces label i by the closest label on the support space Ω:

  Figure 3.8.: Comparison of erosions of black color: (a) original image I, (b) erosion ε black (I; B) using majority-based formulation (3.15), (c) erosion ε black (I; B) using distance-based formulation (3.16). Let l be the width of the red line, the structuring element B is a square whose size is now between l and 2l. (d) iterated black erosion ε 2 black (I; B) using majority-based formulation, (e) erosion ε black (I; 2B) using majority-based formulation.

  Figure 3.9.: (a) Original labeled image I, (b) rst propagation regime (c) second propagation regime.

  Figure 3.11.: n-ary processing of Fig. 3.10 (a). rst row: composed n-ary lter ψ(I; B), second row: iterated composed n-ary lter ψ k (I; B) until convergence, third row: convergence speed w.r.t. to k. Column (a) ψ(I; B) = γ 1 • γ 2 • γ 3 • γ 4 , (b) ψ(I; B) = γ 2 • γ 1 • γ 3 • γ 4 , (c) ψ(I; B) = γ 3 • γ 4 • γ 1 • γ 2 , (d) ψ(I; B) = γ 4 • γ 3 • γ 2 • γ 1 .

  Figure 3.12.: Image size-regularization: (a) original image I, (b) classical order-based ltering, (c) n-ary based ltering, (d) zoom-in (a), (e) zoom-in (b), (f ) zoom-in (c) , (g) residue between (d) and (e), (h) residue between (d) and (f ). See the text for details.

Figure 3 .

 3 Figure 3.13.: Image size-regularization: (a) original image I, (b) classical order-based ltering, (c) n-ary based ltering, (d) zoom-in (a), (e) zoom-in (b), (f ) zoom-in (c), (g) zoom-in (a), (h) zoom-in (b),(i) zoom-in (c). See the text for details.

  us consider the example ofFig.3.15(a), which corresponds to a grey-scale image of a triphasic material. Each phase is characterized by a range of grey level. On a semantic level, each phase has an equivalent level. It is then arbitrary to order them according to their grey level. A fuzzy extension of the n-ary framework is a natural framework to process these data. The restriction to a two-classes mixture is particularly adapted to images originally acquired in grey-scale. The two-classes mixture is built following the fuzzy classication of grey-scale described in the diagram Fig.3.14 for n = 3 phases. Values a, b and c were chosen manually for the experiment but the procedure can easily be automated.

Figure 3 .

 3 Figure 3.14.: Fuzzy classication of the grey levels

55 4 .

 4 Theoretical notions and notations of probability and dierential geometry Résumé Ce court chapitre introduit des notions théoriques et les notations nécéssaires à la compréhension de la deuxième partie de la thèse. Nous avons choisi de ne pas introduire la notion de variété. En eet, mis à part les espaces de rotations, tout les espaces présents dans la suite de l'étude peuvent être décrits en une unique carte: ils sont tous homéomorphes à R n pour un certain n. Il n'y a donc pas besoin d'introduire une dénition abstraite de plan tangent à partir de classes d'équivalences de courbes. Pour plus de détails sur la théorie de la mesure, voir[START_REF] Tao | An Introduction to Measure Theory[END_REF], sur la géometrie Riemannnienne voir[START_REF] Berger | A Panoramic View of Riemannian Geometry[END_REF], et sur l'analyse harmonique dans les espaces symétriques, voir[START_REF] Helgason | Geometric Analysis on Symmetric Spaces[END_REF].

Figure 4

 4 Figure 4.1.: Let X be a random variable valued in [1, 2]. The blue curve is the density of µ X with respect to the Lebesgue measure of R restricted to [1, 2]. The red curve is the density of µ X with respect to the measure µ([a, b]) = log(b) -log(a).

Figure 4

 4 Figure 4.2.: Riemannian space.

  Figure 4.3.: Exponential map.

  Figure 5.1.: Image (a) is a 16 bits gray-scale image. Graph (b) represents the number of pixels associated with each of the 2 16 levels. This is often the quantity used in image processing. Graph (c) is a kernel density estimate of the distribution of the gray values. The kernel bandwidth is determined by a cross-validation technique. Graph (d) is the density computed using the coarea formula. The computation is achieved using an interpolation of the image by pieces of linear functions.

Fig. 5

 5 Fig. 5.2(a) shows a set of ellipses in the projective ab plane of the CIELAB color space. Let us assume that these ellipses represent the local perceptual metric. Let c be a point where the metric has been measured through the ellipse E c . In a neighborhood of c, computing distances using the metric measured at c is a better approximation of the perceptual distance than using the canonical Euclidean distance of the ab plane. At a point p where the metric is originally unknown, a metric interpolated from the neighbor points c i has all the odds of being more relevant than the canonical Euclidean metric of

  dependent of the luminance. Nevertheless, for visualization purpose we choose to abandon the luminance information in order to work with two dimensional data. The MacAdam ellipses were measured at a xed luminance, in the CIE chromaticity diagram. The ellipses are transported in the L = 40 plane of the CIELAB space. Forgetting the luminance coordinate, one obtains then a transport of the MacAdam ellipses in an ab plane. Remind that the proposed framework is independent of the dimension and can be used in three dimensional spaces with standard datasets of ellipsoids.

Fig. 5 .

 5 Fig. 5.3 represents the density of the Riemannian measure with respect to the Lebesgue density of the plane. Recall that the expression of the density is given by det(G) where G is the metric tensor derived from the ellipse.

  Figure 5.2.: (a) MacAdam ellipses transported in the projective ab plane, (b) Delaunay triangulation of the space from center of MacAdam ellipses, (c) interpolated ellipses.

  Assume the term det(G(x)) known for all x ∈ U . According to Chapter 4, dvol dLeb (x) = det(G(x)).

  of the Riemannian measure of M with respect to the Lebesgue measure Leb p after the identication of M and T p M induced by exp p , see Fig.6.1.

Figure 6

 6 Figure 6.1.: Volume change θ x i induced by the exponential map Given K and a positive radius r, the estimator of f proposed by Pelletier (2005) is

Fig. 6 .

 6 Fig. 6.2 sums up the situation. Let ΠA be the restriction of Π A to H A . The density of interest is a product of two factors, i.e.,

Figure 6

 6 Figure 6.2.: The application Π is a Riemannian submersion for the at metric on Gl(n) and the Wasserstien metric on SPD(n). The volume change between the blue and the red area appears in the expression of kernels.

  Fig. 6.3 shows several cases of density change for n = 2.

  Figure 6.3.: For n = 2, (a),(b) and (c) are visualizations of dµ Eucl λ , dµ W ass λ

  Figure 6.4.: (a): Convergence of density estimations using an Euclidean and a Riemannian Kernel, in a Riemannian context, (b): the density is estimated following section 6.2.3 using two standard kernels of the plane, K(||x||) = 3 π (1 -||x|| 2 ) 2 1 ||x||<1 and K(x) = 220 81π (1 -||x|| 3 ) 3 1 ||x||<1 .

  Figure 6.5.: Time lapse sequence of images from the retina: (a) mean image, (b) standard deviation image.

  Figure 6.7.: Densities estimated from the structure tensor eld computed from image (a). Top row, densities computed from the induces measures µ . λ : (b) Euclidean metric of R n × Sym(n), (c) Fisher metric, (d) Wasserstein metric. Bottom row, densities computed using the quotient metrics: (e) Euclidean metric of R n × Sym(n), (f ) Fisher metric (note the change of scale), (f ) Wasserstein metric.

Figure 6

 6 Figure 6.8.: Texture image segmentation: (a) and (b) are respectively a texture image and its associated structure tensor eld. The density of the eigenvalues of the structure tensors is estimated with respect to the two measures discussed in section 6.13 induced by the Wasserstein metric: (c) segmentation with respect to the Lebesgue measure on eigenvalues, and (d) segmentation with respect to µ W ass λ

  Figure 6.9.: DTI segmentation by watershed partitioning of its estimated density (histogram): (a) original SPD(3)-valued image, (b) its map of fractional anisotropy. The density of the largest eigenvalues of each matrix is evaluated against the Lebesgue measure of R n and (c) is a segmentation corresponding to the watershed transform of the complement of the this density.

(

  2013). D endowed with the metric ds is called the Poincaré disk. The application f (z) = z -i z + i , is an isometry between the Poincaré upper half plan and the Poincaré disk. Thus the Poincaré upper half plan and the Poincaré disk are isometric. Fig. 6.10 presents the estimations of the marginal densities coecients µ k . Data used in the experimental tests are radar observations from THALES X-band Radar, recorded during 2014 eld trials campaign at Toulouse Blagnac Airport for European FP7 UFO study (Ultra-Fast wind sensOrs for wake-vortex hazards mitigation). Data are representative of Turbulent atmosphere monitored by radar in rainy conditions. Fig. 6.10 and 6.11 illustrate the density estimation of the six reection coecients on the Poincaré unit disk. For each coecient the dataset is composed of 120 draws. µ 1 , N = 0.18 µ 2 , N = 2.13 µ 3 , N = 4.81 µ 4 , N = 3.83 µ 5 , N = 5.26 µ 6 , N = 6.99

Figure 6 .

 6 Figure 6.10.: Kernel density estimation of the distribution of the 6 rst coecients µ k under rainy conditions, with K(x) = 3 π (1-x 2 ) 2 1 x<1 . For visualization clarity, densities are normalized by N = sup z∈D { f (z)}. Contours correspond to the watershed transform of the complement of estimated densities.

  µ 1 , N = 0.007 µ 2 , N = 1.61 µ 3 , N = 14.86 µ 4 , N = 8.82 µ 5 , N = 17.31 µ 6 , N = 23.29

Figure 6 .

 6 Figure 6.11.: Kernel density estimation of the distribution of the 6 rst coecients µ k in normal conditions, with K(x) = 3 π (1 -x 2 ) 2 1 x<1 . For visualization clarity, densities are normalized by N = sup z∈D { f (z)}. Contours correspond to the watershed transform of the complement of estimated densities.

  les diérents résultats de la thèse et présente les perspectives associées à ces résultats. Le lemme 2.3.1 du chapitre 2 formalise les limites de la qualité des ordres totaux dans les espaces métriques, en terme de morphologie mathématique. A chaque ordre total, il est possible de trouver des images pour lesquelles les érosions et les dilatations sont irrégulières et incohérentes. Il est possible de construire une fonction de coût mesurant la qualité d'un ordre total. La procédure de minimisation de cette fonction de coût proposé dans le chapitre 2 permet de construire un ordre adapté à une image donnée. Les résultats de la procédure de minimisation sont intéressants mais sont encore sans doute sujets à améliorations. Les limites des ordres totaux ont conduit à s'orienter vers un autre type de données : les images à valeur dans un ensemble de labels non ordonnés. Un nouveau cadre morphologique est proposé pour traiter ces images. Une partie importante des propriétés des opérateurs binaire est conservée dans ce nouveau cadre. Les propriétés théoriques étant maintenant en majeure partie déterminées, les travaux futurs se concentreront sur les applications. La formalisation de la notion d'histogramme d'image en terme d'estimation de densité suggère l'utilisation de la formule de la coaire. Cette idée à en fait déjà été explorée dans Duy et al. (2013) dans le cas des images I, R p → R q≤p . Un travail peut-être envisagé sur le cas des images I, R p → R q>p . L'étude des histogrammes d'images couleurs dans le cadre des métriques perceptuelles propose de se baser non pas sur une approximation euclidienne globale de la métrique mais seulement sur une approximation locale. Cela permet d'être plus dèle à la métrique perceptuelle, sans augmenter la complexité des calculs. Ce travail soulève le problème important de l'interpolation d'une métrique riemannienne. Le chapitre 6 passe en revue plusieurs techniques d'estimation de densités de probabilités dans le cas où les données sont des lois Gaussiennes. Parmis les méthodes étudiées, la méthode des noyaux semble la plus adapté à la plupart des situations. Les travaux futurs se concentreront sur les applications pratiques. mathematical morphology led to the following results. Lemma 2.3.1 of chapter 2 formalizes a limitation of the quality of total orders on multidimensional spaces

  faithful to the Riemannian metric without requiring additional computation cost. The analyzed example shows the impact of the local adaptation.Non-parametric probability density estimation has been discussed in Chapter 6 for several standard spaces of Gaussian laws, under dierent metrics. For each case, several density estimation techniques have been discussed and compared. Expressions of kernels for kernel density estimation, and functions involved in orthogonal series density estimation are provided when possible. The most original contributions of this work consists in the expression of kernels for the Wasserstein metric, and in the study of the structure of partial quantities (mean, eigenvalues and rotation) under dierent metrics. Partial quantities play an important role since studied phenomena often depend on only one or two parameters results on mathematical mathematical morphology lead to the following perspectives:

  shift algorithm is very similar to a watershed transform of the estimated density of Morphology, Géométrie et Statistiques en imagerie non-standard.Résumé : Le traitement d'images numériques a suivi l'évolution de l'électronique et de l'informatique. Il est maintenant courant de manipuler des images à valeur non pas dans {0, 1}, mais dans des variétés ou des distributions de probabilités. C'est le cas par exemple des images couleurs où de l'imagerie du tenseur de diffusion (DTI). Chaque type d'image possède ses propres structures algébriques, topologiques et géométriques. Ainsi, les techniques existantes de traitement d'image doivent être adaptés lorsqu'elles sont appliquées à de nouvelles modalités d'imagerie. Lorsque l'on manipule de nouveaux types d'espaces de valeurs, les précédents opérateurs peuvent rarement être utilisés tel quel. Même si les notions sous-jacentes ont encore un sens, un travail doit être mené afin de les exprimer dans le nouveau contexte.Cette thèse est composée de deux parties indépendantes. La première, « Morphologie matématiques pour les images non standards », concerne l'extension de la morphologie mathématique à des cas particuliers où l'espace des valeurs de l'image ne possède pas de structure d'ordre canonique. Le chapitre 2 formalise et démontre le problème de l'irregularité des ordres totaux dans les espaces métriques. Le résultat principal de ce chapitre montre qu'étant donné un ordre total dans un espace vectoriel multidimensionnel, il existe toujours des images à valeur dans cet espace tel que les dilatations et les érosions morphologiques soient irrégulières et incohérentes. Le chapitre 3 est une tentative d'extension de la morphologie mathématique aux images à valeur dans un ensemble de labels non ordonnés.La deuxième partie de la thèse, « Estimation de densités de probabilités dans les espaces de Riemann » concerne l'adaptation des techniques classiques d'estimation de densités non paramétriques à certaines variétés Riemanniennes. Le chapitre 5 est un travail sur les histogrammes d'images couleur dans le cadre de métriques perceptuelles. L'idée principale de ce chapitre consiste à calculer les histogrammes suivant une approximation euclidienne local de la métrique perceptuelle, et non une approximation globale comme dans les espaces perceptuels standards. Le chapitre 6 est une étude sur l'estimation de densité lorsque les données sont des lois Gaussiennes. Différentes techniques y sont analysées. Le résultat principal est l'expression de noyaux pour la métrique de Wasserstein. Mots clés : Morphologie Mathématique, Ordres Totaux, Traitement d'Image, Estimation de Densités, Géometrie Riemannienne, Métrique de Fisher, Métrique de Wasserstein, Histogrammes d'Images. Morphology, Geometry and Statistics in non-standard imaging.

  Select in S 1 the nearest neighbour of Cluster index greater than Cluster index and the nearest neighbour lower than Cluster index . Order elements of Cluster index according to their distance to the two clusters selected during the previous step.

	A function FuncRecursive
	FuncRecursive;
	input : An ordered list of clusters (S 1 , < 1 ) , an Index Cluster Cluster index , a
	depth D
	output: An order on Cluster index given by a list
	if D > threshold then
	SimpleOrder (C luster index )
	else
	S 1 ← Neighbour(S 1 , Cluster index );
	S index ← IndexCutting(Cluster index );
	S 2 ← S 1 ∪ S index \Cluster index ;
	(< 2 ) ← MainOrder(S 1 , Cluster index , S index );
	listorder ← empty list;

  RésuméLa morphologie mathématique pour les images binaires peut être entièrement d'écrite dans le cadre de la théorie des ensembles. Cependant, cette théorie n'est pas susante pour d'écrire la morphologie des images en niveau de gris. Ce type d'images requiert l'introduction de la notion d'ordre. Plus généralement, la morphologie mathématique est maintenant d'écrite à l'aide de la notion de treillis. Depuis quelques décennies, plusieurs

travaux visent à appliquer des opérateurs morphologiques aux images multivariées telles que les images couleurs. Ces travaux ce basent principalement sur la notion d'ordre vectoriels. Cependant, aucune de ces tentatives n'a donné de résulats entièrement satisfaisants.

Au lieu de viser directement le cas multivarié, nous proposons ici une extension de la morphologie binaire à une situation intermédiaire : les images à valeur dans un ensemble ni de labels non ordonnés. Nous proposons ensuite une extension à un cas continu.

Ce chapitre est principalement extrait

de Chevallier et al. (2015a)

.

  .4.2), i.e.,

	G F isher (m,σ) =	1/σ 2 0	0 2/σ 2		(6.35)
	Consider now the application s : (m, σ) → (x, y) = (m,	√	2σ). Solving the equation:
	(ds) t .G F isher (x,y) .ds = G F isher (m,σ)		(6.36)

Acknowledgments

I am convinced I have been very lucky to be supervised by Jesùs Angulo, he has been an ideal director. Thanks to his right balance between guidance and freedom, the three last years have been an exciting scientic experience.

I would like to thank Jasper van de Gronde, for the multiple discussions and constructive suggestions during our respective academic visits. Jasper is always ready to spend time to help his colleagues.

I would like to thank Daniel Bennequin, Jean-Philippe Anker, Frédéric Barbaresco and Ivar Farup for their help on various problems. I also would like to thank the referees for their help in the improvement of the manuscript.

At last, I would like to thank my family, my friends from Malakableau and all the members of the Center for Mathematical Morphology. Thanks to them, the last three years have been particularly memorable.

dened by:

(6.4)

The corrective factor θ x i (x) -1 is necessary since the kernel K originally integrates to one with respect to the Lebesque measure, while we want it now to integrate to one with respect to the Riemannian measure. It can be noted that this estimator is the usual kernel estimator is the case of Euclidean space. Under reasonable assumptions on the true density f , the shape of the kernel does not have a signicant impact on the quality of the estimation in the Euclidean context, see for instance [START_REF] Silverman | Density Estimation for Statistics and Data Analysis[END_REF]. Simulations summarized in Fig. 6.4 experimentally conrms the result in a specic Riemannian situation.

Let r inj = sup p∈M (r inj p ). Let δ = sup K the supremum of the sectional curvature in M. For a denition of the sectional curvature, see. [START_REF] Carmo | Riemannian Geometry[END_REF]. Then, [START_REF] Pelletier | Kernel density estimation on Riemannian manifolds[END_REF] provides the following result. Theorem 6.2.1 [START_REF] Pelletier | Kernel density estimation on Riemannian manifolds[END_REF]) Let x i be an arbitrary point of M. Let µ be a measure on M whose density with respect to the measure vol is

For r < min{ r inj 2 , π 4 √ δ } when δ > 0 or r < r inj 2 when δ ≤ 0, x i is an intrinsic mean of µ.

Theorem. 6.2.1 ensures the kernel put on x i is centered on x i .

Despite convergence rates provided by [START_REF] Pelletier | Kernel density estimation on Riemannian manifolds[END_REF] are obtained for compact manifolds without boundaries, they remain valid for non-compact manifolds when the support of X is compact. Indeed the double manifold, see [START_REF] Munkres | Elementary Dierential Topology[END_REF], enable to see a compact manifold with boundaries as a submanifold of a compact manifold without boundaries.

Assume f two-times dierentiable with bounded second covariant derivative. See [START_REF] Parker | Geometry Primer[END_REF] for a denition of covariante derivative. [START_REF] Pelletier | Kernel density estimation on Riemannian manifolds[END_REF] also provides the following convergence rates.

Theorem 6.2.2 [START_REF] Pelletier | Kernel density estimation on Riemannian manifolds[END_REF]) For r satisfying conditions of theorem. 6.2.1, there exist a constant C f such that:

For the expression of the constants see details in [START_REF] Pelletier | Kernel density estimation on Riemannian manifolds[END_REF].

Given a reference point p ref ∈ M, let

Theorem 6.2.5 [START_REF] Hendriks | Nonparametric Estimation of a Probability Density on a Riemannian Manifold Using Fourier Expansions[END_REF]) Suppose f is s times dierentiable with square integrable derivatives. Let T 0 > 0. There are constants A (depending on M) and B (depending

and for a suitable decrease of T ,

).

Thus, the convergence rate for the orthogonal series is similar to the rate of the kernel estimation.

The Fourier series on a parallelepiped of R n are a particular case of decomposition on eigenfunctions of the Laplacian.

When the space is not compact, the estimator takes an integral form. In the case of a continuous spectrum,

for some uncountable set S. Then, in this uncountable case, one has

where d * s is the spectral measure on S and in the case of symmetric spaces is the Harish-Chandra c-function, see [START_REF] Helgason | Dierential Geometry, Lie Groups, and Symmetric Spaces[END_REF].

The explicit expressions of F T and of the estimators are provided when necessary in the following sections. For the space of positive denite matrices equipped with the metric invariant by the action of Gl(n) on SPD(n),(G, M ) → GM G

H(∆f

Then, ∆ σ/2 is dened as the operator verifying:

Note that H is a group. The quotient SO(n)/H = {RH, R ∈ SO(n)} can be identied to the fundamental domain D ⊂ SO(n) where for R ∈ D, the maximal element with the smallest index element of each column 1 ≤ j < n is positive. In other words

Then, parametrization ψ :

is an homeomorphism.

Let α be the change of parametrization α = φ -1 • ψ:

On any locally compact group there exists, up to a scaling factor, a unique measure invariant under the group law. Such a measure is called a Haar measure. The Lebesgue measure is a Haar measure. Since the set of rotations is a compact group, there exist a

Haar measure Haar SO(n) on rotations. Let µ φ be the Lebesgue measure on

Let µ ψ be the product measure between the Lebesgue measure of R n × R n + and the Haar measure of SO(n), i.e.,

We are now interested in the local volume change induced by the application α, that is to say the ratio between the volume of an innitesimal volume element in

. By denition of the Haar measure, µ ψ (A) = µ ψ (R.A). Thus the ratio between the volumes of innitesimal volume elements are independent of the rotation.

denotes the density of µ φ with respect to α * (µ ψ) at (m, Σ). For Σ = RD λ R t with D λ ∈ E, R ∈ D, we abuse notation and write

(m, D λ , I)

Separation of mean, eigenvalues, and rotation

This section shows that in many classical cases, measures on multivariate Gaussian laws can be decomposed as a product between a measure on the mean space, a measure on the eigenvalue space, and a measure on the rotation space. Observe that rotations naturally

Theorem 6.4.1 Let µ be a measure on G(n), invariant under the action of rotations. Measure µ is equivalent to a product measure,

where

Proof This is a direct consequence of Theorem 8.2, page 155, of [START_REF] Wijsman | Invariant measures on groups and their use in statistics[END_REF].

In most situations, measures µ on G(n) present the desired invariances, and µ(G(n) \ G(n)) = 0. Writing µ as a product measure enables to dene densities for the mean, the eigenvalues and the rotations separately, which is not always possible, see Section 6.13. This is in particular the case for the measures induced by the Fisher metric, the Wasserstein metric, and the Lebesgue measure of R n × SPD(n).

The Fisher metric

Unlike multivariate centered and univariate Gaussian laws, the space of multivariate Gaussian laws under the Fisher metric is not a symmetric space. As in the large majority of the Riemannian manifolds, there is no explicit expression of distances. We only have access to the expression of the metric. An attempt was made by authors of [START_REF] Lovri¢ | Multivariate Normal Distributions Parametrized as a Riemannian Symmetric Space[END_REF] to symmetrize the space.

Let u and v be vectors of R n Sym(n). Let index m and Σ denote the components associated to R n and Sym(n) respectively.

The scalar product of the Fisher metric at (m, Σ) can be rewritten, see [START_REF] Skovgaard | A Riemannian Geometry of the Multivariate Normal Model[END_REF],

be the matrix of the scalar product in the canonical basis of R n Sym(n). Let vol F isher be the measure associated with the Riemannian metric. The density of the Fisher volume measure with respect to the Lebesgue measure of R n × SPD(n) is given by

The invariance of the Fisher metric under the action of rotations implies the invariance of the measure, i.e.,

Thus the volume change only need to be computed for the case diagonal matrices. Let D ∈ SPD(n) be a diagonal matrix of eigenvalues λ i . Using Expression (6.11), we can now For s ∈ C n and K ∈ SO(n), the Helgason functions (plane waves) have the following expression:

Let C ∞ c be the set of innitely dierentiable functions with compact support. The FourierHelgason transform of f ∈ C ∞ c is given by:

Σ)dvol F isher , (6.27)

where x denotes the complex conjugate of x.

where Γ is the gamma function. For s

,

and

Using this notation, the measure µ s on C n is dened as Terras (1988) and Peter T. Kim and Donald St. P. Richards (2011) provided the inversion formula,

where dk stands for the normalized factorization of the Haar measure of SO(n) on SO(n)/H.

Following the discussion of section 6.4.1, the expression can be rewritten as

where Haar SO(n) is taken such that Haar SO(n) (D) = 1.

Density estimation on multivariate centred Gaussian laws

Let X be a random variable in Ω → G m=0 (n). Notations of this section are similar to those of section 6.2. For the reasons explained in the previous section, we are here only interested in the Fisher metric. The case of the Wasserstein metric is the same studied in section 6.5.2.

Histograms

The problem has already been partially addressed in section 6.7.1. Recent works on Voronoï diagrams in the hyperbolic space Nielsen and Nock (2010) might provide alternative solutions through the construction of adaptive bins.

Kernel

The explicit expression of the change of volume θ p (q) induced by the exponential map enables us to use the kernel density estimation. Given a scaling parameter r > 0, the estimation becomes: (6.46) where θ z i (z) is given in Eq. ( 6.42). Then, we have the close expression: (6.47) where here d F isher (z,

Orthogonal series

The symmetry of the space provides a Helgason transform and its inverse. The empirical transform is given by:

where e s,ku is the plane wave dened in 6.43. For a cut-o parameter T > 0, the associated density estimator is then given by:

and, consequently the estimation do not privilege specic directions or locations.

The main drawback of this method is its computational cost. Indeed the situation is similar to that considered in section 6.7.3.

Partial quantities: mean, eigenvalues and rotation

This section addresses the study of the standard partial quantities that are the mean, the eigenvalues and the rotation. On the one hand, depending on the application, the interesting information is sometimes carried by only one or two partial quantities. On the 6.13.2. Quotient metric on partial quantities

The space associated to each partial quantity can be seen as a quotient of G(n) by the rest of the partial quantities. The quotient G(n)/(R n + ×SO(n)) denotes the identication elements with all those which dier only in the eigenvalues or the rotation. Thus G(n)/(R n + ×SO(n)) can be interpreted as the space of means. In several cases, the Riemannian metric on G(n) induces a canonical metric on such quotient space. This metric induces a measure on the partial quantity, which enable the estimation of density of the associated random variable X (.) .

Mean

For the Wasserstein metric and the Euclidean metric of R n × Sym(n), the quotient space R n ∼ G(n)/(R n + × SO(n)) inherits naturally of the canonical Euclidean metric. Indeed, each couple (λ, R) induces the same metric on m. The density estimation can thus be achieved using standard techniques. Note that the induced Lebesgue measure is similar to µ mean .

Recall the expression of the scalar product of the Fisher metric (6.11):

Each covariance matrix Σ induces a metric on the space of means. However, the metric varies when Σ varies. Thus the Fisher metric does not induces a canonical metric on

Eigenvalues

For λ ∈ R n , let D λ be the associated diagonal matrix. Fixing a mean m and a rotation R induces a metric on eigenvalues:

(6.51)

The group of rotations naturally acts on Gaussian laws by R a . Lef t N m,λ,R b = N m,λ,RaR b . Furthermore for the Euclidean metric of R n × Sym(n), the Fisher and the Wasserstein metric, we have that for all R ∈ SO(n):

Thus since d m,R is independent of m and R, the three metrics induce a canonical metric on the quotient R n + ∼ G(n)/(R n × SO(n)).

On the one hand according to section 6.8, the Wasserstein metric on eigenvalues at xed rotation is the canonical Euclidean metric after the change of coordinates σ = √ λ. On the other hand, according to section 6.9 the Fisher metric on eigenvalues at xed mean and rotation is the canonical Euclidean metric after the change of coordinates γ = log(λ).

Finally, it is easy to see that the quotient metric induced by the Euclidean metric of R n × Sym(n), is the Euclidean metric on R + . Note that the associated measures diers from µ λ .

Conclusions and perspectives

The whole theory and corresponding practical solutions to density estimation have been analyzed for standard spaces underlying Gaussian laws.

Under each space and each metric, four types of density estimation are considered.

The rst one consists in classical estimations in the Euclidean context of the parameters, multiplied by a density ratio to obtain a density with respect to the desired measure. This type of density estimation is possible in each studied space and metric, but does not take into account the geometry of the space: if two metrics induce the same volume measure then the estimated density does not depend on the choice of one the metrics. The three other types of estimation, histograms, kernels and orthogonal series rely more directly on the geometry. On the one hand, despite their low computational cost, nding good tiling make the use of histograms dicult. On the other hand, nding an orthogonal basis of functions adapted to orthogonal series density estimation is also a dicult problem. In the few situations where the eigenfunctions of the Laplacian operator are known, the estimator presents a high computational complexity.

The kernel density estimation seems to be the most adapted in most cases. It presents a reasonable computational complexity, and expressions of kernels are ready to use in each studied situations, except the space of multivariate Gaussian laws under the Fisher metric.

To our knowledge, the most original contribution of this work are the expression of kernels for the case of the Wasserstein metric, and the study of partial quantities. The mean shift algorithm is a standard approach that search maxima of a density associated with a point cloud and then segment it into modes or classes. This algorithm, proposed by [START_REF] Fukunaga | The estimation of the gradient of a density function, with applications in pattern recognition[END_REF] has been widely studied in R n . It is mainly used in 7.2. Perspectives the point cloud [START_REF] Paris | A Topological Approach to Hierarchical Segmentation using Mean Shift[END_REF]. The theoretical properties of the algorithm on spaces of Gaussian laws and its applications form an interesting subject of future studies.

Density estimation is also an important tool in classication. This is the subject of recent studies on symmetric positive denite matrices using parametric density under the ane-invariant metric. Recall that the ane-invariant metric is similar to the Fisher metric on multivariate centered Gaussians. In the recent work [START_REF] Said | Riemannian Gaussian Distributions on the Space of Symmetric Positive Denite Matrices[END_REF], it is proposed to model the distribution of each classes by a parametric model and to take classication decision according to the estimated probabilities. The expression of kernels for the Fisher and the Wasserstein metric enables to conduct similar studies using non-parametric density estimation.