Je Remercie

Isabelle Bloch D'avoir

Michel Bilodeau

Serge Beucher

Etienne Decencière

Petr Dokladal

Dominique Jeulin

Beatriz Marcotegui

Fernand Meyer

François Willot

Bruno Figliuzzi

Santiago Velasco-Forero

Serge Koudoro

Matthieu Faessel

Un

Aurélie, Bassam, Benjamin, Daniele, Elise, Emmanuel, Enguerrand, Haisheng, Joris Ça Je Remercie Amin

Jean Baptiste

Lydia Luc

Vaïa Vous Avez

Giancarlo Guillaume

Jeremy Laura

Marco Laurent

Keywords: Image processing, Machine learning, Kernel methods, Mathematical morphology, Principal Component Analysis, Support Vector Machine, Deep learning, Scattering transform, Kriging Mots Clés Image processing, Machine learning, Kernel methods, Mathematical morphology, Principal Component Analysis, Support Vector Machine, Deep learning, Scattering transform, Kriging

J'aimerais tout d'abord remercier mon directeur de thèse, Jesus Angulo, pour son encadrement hors pair. Ainsi que pour la confiance qu'il m'a accordée en acceptant de m'encadrer, et pour toutes les heures qu'il a consacrées à me diriger, à m'aider et à me conseiller durant mes recherches. Sa grande disponibilité et ainsi que nos nombreux échanges m'ont permis de mieux définir mon projet. Je souhaite également remercier Maxime Monreaud et Loïc Sorbier, mes encadrants à IFP Energie Nouvelles, pour leur aide et leurs conseils tout au long de ma thèse, ainsi que leur implication tout au long de mon projet qui a permis de rendre possible beaucoup de nos idées.

Ce travail n'aurait pas été possible sans les Mines ParisTech, IFP Energie Nouvelles ainsi que de l'association des anciens des Mines ParisTech, qui grâce à leur soutien financier m'ont permis de me consacrer sereinement à mes travaux de recherches.

Je remercie également les membres de mon jury de thèse pour avoir formulé des critiques constructives et pour leurs conseils.

Résumé

Cette thèse porte sur la statistique spatiale multivariée et l'apprentissage appliqués aux images hyperspectrales et multimodales. Les thèmes suivants sont abordés :

Fusion d'images :

Le microscope électronique à balayage (MEB) permet d'acquérir des images à partir d'un échantillon donné en utilisant différentes modalités. Le but de ces études est d'analyser l'intérêt de la fusion de l'information pour améliorer les images acquises par MEB. Nous avons mis en oeuvre différentes techniques de fusion de l'information des images, basées en particulier sur la théorie de la régression spatiale. Ces solutions ont été testées sur quelques jeux de données réelles et simulées.

Classification spatiale des pixels d'images multivariées :

Nous avons proposé une nouvelle approche pour la classification de pixels d'images multi/hyper-spectrales. Le but de cette technique est de représenter et de décrire de façon efficace les caractéristiques spatiales / spectrales de ces images. Ces descripteurs multi-échelle profond visent à représenter le contenu de l'image tout en tenant compte des invariances liées à la texture et à ses transformations géométriques.

Réduction spatiale de dimensionnalité :

Nous proposons une technique pour extraire l'espace des fonctions en utilisant l'analyse en composante morphologiques. Ainsi, pour ajouter de l'information spatiale et structurelle, nous avons utilisé les opérateurs de morphologie mathématique.

Table of Contents

I Introduction

II Feature representation and classification for hyperspectral images

Multivariate images

A grey scale image can be seen as a matrix in which each pixel is an entry of the matrix. Such an image brings information on the spatial structures present in it; a natural evolution was to add the color information. On color images each pixel is a vector of size 3. Thanks to these images we have both spatial and spectral information. However the spectral information is quite limited since we just have access to the colors visible by human eye. This can be sufficient in some cases, but if we want to be more precise we might need more spectral information. That consideration gave birth to hyperspectral images. On each pixel of these images we have a vector which is a signal that represents the reflected light of the area under study. Color or hyperspectral images are examples of multivariate images. Multivariate images are images where each pixel is described by a vector of size superior to 1, contrary to traditional univariate (or scalar) images. These types of images are becoming more important because of the sensors abilities and can arise from a huge variety of sources. Multivariate images appear in remote sensing, astronomy, molecular imaging, etc. The image analysis of these sets must be different, and adapted to the nature of the data itself. In this thesis we will work with two kind of multivariate images. On the one hand, on multispectral and hyperspectral images in Part 1. Typical examples of these images are obtained by remote sensing devices, and our objective is mainly to classify each pixel of them, using the spectral and spatial information. The second kind of data we are interested in is the Energy Dispersive Spectrometry (EDX), that are used in physicochemical characterization of materials. The work on these images are the fruit of a collaboration between MINES ParisTech and IFP Energies Nouvelles.

Curse of dimensionality

The curse of dimensionality is a term invented by Richard Bellman [13] which designates the phenomena that occur when data are studied in a space of high dimension. Suppose that the data can be represented by a vector v ∈ [0, 1] D , where D is the dimension of the data space. Consider first a simple case where D = 1, in this case it is very easy to calculate that with 100 points, you can get an interval between the points which is about 10 -2 . However, if we consider a space of dimension D = 10 and if you want the items to be separated by a ball of radius 10 -2 , we easily see that it takes 10 20 points. Now, if we consider a space of dimension D, we see that we need 10 2D points which can be astronomical. This shows that to fight the curse of dimensionality we must increase the number of data. However most of the time the number of data is limited. Let us consider that we have 900 vectors v i i ∈ [1, 900] of dimension 2, see example in Figure 1(a). These vectors are the result of a Gaussian mixture model, where three Gaussians have been considered. To study these vectors, we plot in Figure 1(b) the distance matrix d, which is a square matrix of size 900 such that d i,j = v iv j 2 . Then, using the k-means clustering algorithm [START_REF] Forgy | Cluster analysis of multivariate data: efficiency versus interpretability of classifications[END_REF], we can naturally separate the data into three clusters, which is illustrated in Figure 1(c). When the dimension of the data is equal to 2, it seems easy to analyze the data. However if the dimension of the data is increased, as in Figure 2(a)(b), we can see that the distance matrix d is not relevant any more for D = 100 and so the k-means is not able to find the three clusters.

v i ∈ [0, 1] D , i = 1, • • • , 500
, where D is the dimension of the feature space. We plot R = max (i,j) { v i -v j 2 } min (i,j) { v i -v j 2 } , which represents the power of discrimination of the distance. To have consistent results, 100 Monte Carlo simulations have been generated, and the black curve represents the mean result, and the dash-dotted lines represent the mean +/-the standard deviation.

Moreover, it can be shown (see [START_REF] Friedman | The elements of statistical learning[END_REF]) that because of the dimension D of the space, the ratio between the maximum Euclidean distance and the minimum Euclidean distance, i.e., R = max (i,j) { v iv j 2 } min (i,j) v iv j 2 } , of the data tends to 1 when D increases. Hence, the euclidean distance has a very poor discriminating power with high dimensional data. We represented this lost of discriminant power of the L 2 norm in Figure 3. Indeed, in the case of ultraspectral images, D is usually about one or two thousand, which shows how sparse the sampled manifold is (or how empty the whole space is), and how crucial is the question about similarity between spectra. Therefore, due to the high dimensionality of the space, it is critical to learn a good representation in order to fight against the curse of dimensionality.

There is the range of approaches belonging to the dimensionality reduction paradigm, so one could work on a smaller dimension space and thus circumvents the curse of dimensionality. It is also possible to work on techniques able to learn the feature space such as kernel learning methods, or the Convolutional Neural Network (CNN) deep learning approaches.

Multimodal images

Sometimes the acquisition of samples can be done using simultaneously different techniques, which leads to the notion of multimodal image acquisition. These different modalities might bring different information sources that could be combined. This is the case for instance in satellite remote sensing (panchromatic images + ultra/hyper-spectral images), in biomedical microscopy (multispectral images + Raman spectroscopy image + quantitative phase) or in scanning electron microscopy (secondary electrons/backscattered electrons + EDX). The goal of multimodal information combination can be approached by different information fusion techniques. This kind of data brings new issues; firstly, they must be perfectly registered or in other words, these various images should be spatially consistent: pixel locations from one modality should be related to the other one. In this thesis, we consider that the data are perfectly registered. Secondly the information fusion should be compatible in terms of their ranges. Finally, some modalities might see objects absent in others.

Context of the study on multimodal SEM

IFP Energies Nouvelles is a public-sector research institute, active in the fields of energy, transport and environment. Among the many areas of expertise, IFP Energies nouvelles contributes to mastering and developing microscopic analysis techniques. Energy Dispersive Spectroscopy (EDX) images are more and more important in microanalysis. With this type of picture, it is possible to see in each pixel the composition of the pixel, and thus theoretically to have the composition of the analyzed sample. This type of image could therefore greatly improve the microscopic analysis of a sample by adding new information. It should be noted, that the resolution of these are often poorer than other images from Scanning Electron Microscopy (SEM), and do not take into account the sample 3D structure. So in order to be able to help IFP experimenters to characterize these images we worked, in collaboration with an IFP Energies Nouvelles team, on techniques able to merge the different modalities of the SEM images. The literature on multimodal SEM fusion is, to the best of our knowledge, very limited. Hence, we took inspiration for our techniques from hyperspectral remote sensing images pansharpening techniques.

Spatial representation

We also worked on the classification problem of remote sensing images, which is an important issue on remote sensing.

Multivariate imaging allows us to obtain typically a vector describing each pixel of the image, and thus with this quantity of information acquired, we can in theory determine more easily the class of each pixel. However, as we have discussed above, a major problem in classification is the curse of dimensionality. It is therefore necessary to use techniques of statistical learning in order to extract useful information. At the beginning, linear machine learning technique were the first kinds of techniques used in remote sensing [START_REF]Hyperspectral Imaging techniques for Spectral Detection and Classification[END_REF][START_REF] Qian Du | Article: Independent-component analysis for hyperspectral remote sensing imagery classification[END_REF]. Unfortunately the nonlinear nature of the data is then completely neglected, so many approaches based on manifold learning methods [START_REF] Tenenbaum | A global geometric framework for nonlinear dimensionality reduction[END_REF][START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF][START_REF] Belkin | Laplacian eigenmaps and spectral techniques for embedding and clustering[END_REF][START_REF] Donoho | Hessian eigenmaps: Locally linear embedding techniques for high-dimensional data[END_REF][START_REF] Coifman | Diffusion maps[END_REF][START_REF] Schölkopf | Kernel principal component analysis[END_REF] were developed. However, all these techniques often neglect the interactions between a pixel and its neighbours, which is an essential concept in image processing. Also, these techniques have often been invented for data that are not images and were created to address the problems of "big data" that are nonlinear. Therefore gradually a large number of works have attempted to use spatial descriptors giving birth to spatial machine learning, more naturally adapted to image classification. Among the descriptors based on a combination of machine learning and image processing, we note the use of mathematical morphology techniques on hyperspectral images to extract interesting features [START_REF] Ouzounis | Differential area profiles: Decomposition properties and efficient computation[END_REF][START_REF] Fauvel | Spectral and spatial classification of hyperspectral data using svms and morphological profiles[END_REF][START_REF] Dalla Mura | Morphological attribute profiles for the analysis of very high resolution images[END_REF]. That brings now one of the best classification results on remote sensing classification. This justifies our choice to work on spatial machine learning.

Main contributions and Thesis outline

The research work presented in this thesis investigates two main topics. The first subject consists in learning appropriate representations for multivariate images. We worked on both dimensionality reduction and classification techniques. The second theme is more domain applied and deals with building efficient representations for the fusion of multimodal SEM images.

This thesis is composed of a series of chapters, most of them are based on the publications done during the three years, and organised as follows:

Part 1: On this part we focus on spatial dimensionality reduction and supervised classification with examples from remote sensing images. First, in Chapter 2 we present techniques able to reduce the dimensionality of the spectra while considering the spatial information. These techniques lean on morphological techniques and are able to learn a manifold that brings spatial and spectral information. In Chapter 3, in order to improve classification of hyperspectral pixels, we propose a deep texture descriptor model able to characterize the hyperspectral in homogenous local areas of the image. This descriptor is based on new innovative techniques, namely the scattering transform [START_REF] Bruna | Invariant scattering convolution networks[END_REF][START_REF] Mallat | Group invariant scattering[END_REF] and the kernel mean map [START_REF] Smola | A Hilbert space embedding for distributions[END_REF].

Part 2: The second part focusses on multimodal SEM image fusion. Chapter 4 provides a background on the EDX multimodal SEM images. Then, in Chapter 5, we present first some state-of-the-art techniques of pansharpening, and second, after applying these techniques and based on the results, we introduce two techniques based on the bilateral filter and on morphological down/upsampling theory. Finally, we propose in Chapter 6 a more theoretical study of image fusion based on kriging, also linked with Gaussian process regression.

Part II

Feature representation and classification for hyperspectral images

Chapter 1

State-of-the-Art on Statistical Learning Applied to Hyperspectral Images

Abstract

In machine learning, high-dimensional data correspond to points lying in spaces of dimension typically higher than 10 and are difficult to represent and so to interpret. Due to the curse of dimensionality, feature extraction and classification in such spaces is challenging. A typical approach consists in "simplifying the data", usually by projecting the data onto a space of low dimension. This kind of techniques is commonly followed in hyperspectral image processing. These images come from hyperspectral sensors and can be seen as a three-dimensional hyperspectral data cube where each image represents a spectral band. In general the number of bands is around ten for multispectral images, about one hundred for hyperspectral images, and of thousands for ultraspectral images. So it is necessary to use statistical learning techniques to have a better understanding of this data. In this chapter we will focus on a review of the statistical learning techniques that have been developed, on the one hand, to reduce the dimension of our data, and on the other hand, to classify pixels from these images. Many different dimensionality reduction techniques have been proposed, we focus here on the Principal Component Analysis (PCA), Multi-Dimensional Scaling (MDS), Isomap, Kernel-PCA, Diffusion Maps, Locally linear Embedding (LLE), Laplacian Eigenmaps and Hessian LLE. We also remind some dictionary learning algorithms. Finally we explain linear regression, linear classification, Support Vector Machine (SVM), and Convolutional Neural Network (CNN) for classification.

Résumé

L'étude des données en grande dimension (supérieur à 10) est difficile en raison du fléau de la dimensionnalité. Ainsi l'extraction d'information utile est dure. Une approche typique consiste à simplifier les données en projetant les données sur un espace de plus petite dimension, on parle alors de réduction de la dimension.

Introduction

Conventionally, in remote sensing, and in other image processing domains, we use instruments whose measurements are in optical, infrared and radar electromagnetic spectrum areas. Hyperspectral imaging is an evolution of optical imaging to reconstruct the spectral profiles of objects imaged by the acquisition of several tens or several hundred narrow spectral bands, that totally covers the whole of the optical spectral range, from purple to infrared (in most of the case but we can choose another part of the spectrum). This type of imaging, because of its complexity and its data size, is still limited to mainly experimental applications. However, progress in technology and information begin to turn it more operational. That is why several space agencies develop hyperspectral sensors. Among many other initiatives, we can highlight the current German and Italian programs (EnMAP and PRISMA), the American missions in progress Hyperion and the new program HyspIRI (coming soon), the SPECTRA program for ESA. But these kinds of images provide rich information that there can also be used in medical imaging, in physicochemical imaging, etc. Let us formalize the models for hyperspectral image data representation. We review dimensionality reduction and classification techniques relevant for hyperspectral image analysis.

Hyperspectral images

In 1666, Isaac Newton showed that light can be decomposed into a spectrum of light rays that can be represented in a graph linking the intensity and the different wavelengths present in a light beam. The objects that we see have colors because they absorb certain wavelengths and reflect others. This reflected light represents what is seen by our eyes. For example: a red apple reflects mainly red light and grass mainly reflects green light.

In the context of imaging, it is the spectrum of reflected light that interests us. However, in traditional imaging we take the wave included in the visible range. Visible light is the part of the electromagnetic radiation that our eyes capture. In fact, the human eye perceives the waves electromagnetic whose wavelengths lie between about 380 and 780 nanometres. However, in the case where the incident CHAPTER 1. STATE-OF-THE-ART ON STATISTICAL LEARNING APPLIED TO HYPERSPECTRAL IMAGES 13 light comes from the sun, by limiting ourselves to the visible spectrum we do not have access to all available information. Each spectrum has a different name and a different function, such as the spectral range of X-rays used in medical radiology. Thus, depending on its composition, an object absorbs such and such wavelengths and returns a spectrum. This spectrum of reflected light is unique and depends on the composition of the object, and is called spectral signature. Knowledge of the spectral signature of a body can help us to know if it is present or not in a region. It is therefore possible to have access to the material composition.

A hyperspectral image is therefore a simultaneous acquisition of images of many narrow, contiguous spectral bands. This technique enables us to obtain a spectrum for each pixel of the image.

Representation of hyperspectral data Physical description

In this section we describe several physical issues that can deteriorate the acquisition of the data, most of them are mainly related to remote sensing hyperspectral imaging, but some issues are more general.

Atmospheric effects During its path between the object and the satellite, each reflected ray passes through the atmosphere, which changes the reflected spectrum. Indeed, the atmosphere is anisotropic, so there may be interactions between particles within it, through absorption and diffusion phenomena. Therefore, the spectrum recovered from the satellite is not the spectrum of the real object. This is why it is important to consider this aspect, which is one of the main obstacles to hyperspectral satellite imagery. It can be seen by comparing the spectra of Figure 1.2 that the spectrum is actually modified.

Topographic effects Topographic effects are due to the shape of objects. For example in Figure 1.3, we see that the rays A, B and C come from different reliefs and are not reflected in the same place. Thus, it is easy to understand that when a sensor receives a signal from a non-planar surface, a part of the signal received by the sensor does not come from the analyzed area.

1.2. HYPERSPECTRAL IMAGES Figure 1.2: Reflected spectra in the left with no atmosphere, in the right in presence of atmosphere [START_REF] Exelis | Analyse spectrale avec ENVI[END_REF].

Figure 1.3: Reflected spectra from a non-planar object [START_REF] Exelis | Analyse spectrale avec ENVI[END_REF].

Spectral Mixture effects When a photo is taken, the picture is divided into pixels. So that each one has a different color. Even if a pixel is composed of different elements, it will have only one color. However, in hyperspectral images, when a pixel (which is a vector) is composed of different materials we have a spectral mixture. In fact, the spectrum of each pixel is a function of the spectra of different materials integrated in the surface of the detector. This type of problem can make impossible to determine the composition of a pixel.

There are two main models of spectral mixture, Linear Mixing: in this case we consider that the materials of a pixel are optically separated, thus the spectrum of the pixel is just the sum of the different spectra multiplied by a coefficient related to the proportion materials. In the following we refer to this coefficient as the abundance ratio or abundances.

Nonlinear Mixing : in this case we are dealing with a mixture of different materials that creates dispersion effects. The signal thus obtained is a non linear combination of the different spectral reflections. There are different non linear mixture models according to the non linear aspect we want to represent.

Figure 1.4: The two models of spectral mixing [START_REF] Eches | Méthodes Bayésiennes pour le démélange d'images hyperspectrales[END_REF].

CHAPTER 1. STATE-OF-THE-ART ON STATISTICAL LEARNING APPLIED TO HYPERSPECTRAL IMAGES 15

Mathematical description of linear mixing model

Despite the above description, which shows that in a realistic situation, the received signal by the hyperspectral sensors is nonlinear, we will consider a linear model, as it turns out to be frequently used. This assumption may seems very simplistic, but it is one of the most used models and gives consistent results. Finally, this model also has the advantage of being simple from a mathematical point of view.

The linear mixing model assumes that the D-dimensional spectrum pixel vector v is a linear combination of the R pure spectra, called endmembers, tainted by an additive noise N .

v = R k=1 α x,k φ k + N , (1.1)
where φ k ∈ R D for k ∈ [1, R] represent the R endmembers, and α x,k represents the abundance ratio of the endmember φ k over v = f (x) which are constrained:

• all the abundance ratios are positive, i.e., α x,k ≥ 0, ∀k, and

• they belong to the simplex, i.e., R k=1 α x,k = 1 ∀x ∈ E.

State-of-the-art on dimensionality reduction 1.3.1 Notation

We introduce here the notation used in the rest of the chapter. Let E be a subset of the discrete space Z 2 , which represents the support space of a 2D image and F ⊆ R D be a set of pixels values in dimension D. Hence, it is assumed in our case that the value of a pixel x ∈ E is represented by a vector v ∈ F of dimension D. This vector v represents the spectrum at the position x. We write f : E → F for the function associated to the image, so we have f (x) = v. Let F be a set of n spectra such that F = {v i } n i=1 ∈ R D and let us assume that these points lie on a smooth manifold F of intrinsic dimension d with d D. We can consider that F is an hyperspectral image, by doing that, we did not take into account the position of the pixel.

Additionally, we denote the manifolds with a calligraphic upper-case letter (I, S, . . .), and D refers to the dictionary in dictionary learning techniques.

From PCA to manifold learning

The goal of manifold learning techniques is to find a smooth mapping Ψ : F -→ F where F ⊂ R d represents some parameter domain. Mapping Ψ must be found while preserving some characteristic properties. Manifold learning techniques can be mainly divided into two families: i) linear and ii) nonlinear methods. Moreover, the nonlinear techniques can also be divided into global methods that uses only global information from the set of points, and the approaches using local information. We present some of the most popular manifold learning techniques. In particular, we provide more details for PCA and Kernel-PCA, together with MDS.

STATE-OF-THE-ART ON DIMENSIONALITY REDUCTION

Principal Component Analysis (PCA) -A linear method

Principal Component Analysis (PCA) has several names, so it is also known as the Karhunen-Loève transform, the Hotelling transform, the Singular Value Decomposition (SVD) transform, etc. Here we will call it PCA.

We start with a set of n points F = {v i } n i=1 ∈ R D . The PCA goal is to reduce the dimension of this vector space finding the basis that captures most of the variance of data set thanks to a projection on the principal component space, namely

F = {v i } n i=1 -→ F = {v i } n i=1 (1.2) with v i ∈ R d , where d D.
A more geometrical point of view of the PCA is to consider that the PCA goal is to find a space of projection such that the mean squared distance between the original vectors and their projection is as small as possible, this is equivalent to finding the projection that maximizes the variance.

Let us call w j ∈ R D the j principal component. The aim of PCA is to find the set of vectors {w j , 1 ≤ j ≤ D} such as: arg min

w j n -1 n i=1 v i -< v i , w j > w j 2 , ∀1 ≤ j ≤ D. (1.3)
Developing now the distance we have:

v i -< v i , w j > w j 2 = v i 2 -2 < v i , w j > 2 + < v i , w j > 2 w j 2 ,
then adding the additional constraint that w j 2 = 1, replacing in (1.3) and keeping only terms that depend on w j , we have the following new objective function:

arg max w j , w j 2 =1 n -1 n i=1 < v i , w j > 2 , ∀1 ≤ j ≤ D. (1.4) Since var(< v i , w j >) = n -1 n i=1 (< v i , w j >) 2 -(n -1 n i (< v i , w j >)) 2 ,
if we consider that the data F has been column-centered, which means that n i=1 v i = 0, then

var(< v i , w j >) = n -1 n i=1 (< v i , w j >) 2 .
Thus we can see that the goal of the PCA is to find principal components that maximize the variance. The problem can be rewritten in a matrix way using:

n -1 n i=1 < v i , w j > 2 = n -1 (F w j ) T (F w j ) = w T j (n -1 (F T F ))w j = w T j V w j , where V = n -1 (F T F ), V ∈ M D,D (
R), is the covariance matrix of F . Hence we should optimize:

arg max w j , w j 2 =1 w T j V w j , ∀1 ≤ j ≤ D.
(1.5) CHAPTER 1. STATE-OF-THE-ART ON STATISTICAL LEARNING APPLIED TO HYPERSPECTRAL IMAGES 17 Thanks to Lagrange multiplier theorem we can rewrite the objective function as:

L(w j , λ) = w T j V w j -λ(w T j w j -1), (1.6)
where λ ∈ R. Since we want to maximize this function, we have to derive it and equal it to zero:

∂L ∂w j (w j , λ) = 2V w j -2λw j = 0.
So, we finally obtain as solution V w j = λw j .

(1.7)

Thus, the principal component w j that satisfies the objective function is an eigenvector of the covariance matrix V , and the one maximizing L(w j , λ) is the one with the larger eigenvalue. Then we can have all the w j by computing the SVD of V . There are different approaches to choose the reduced dimension d, that is the principal component to be kept. The underlying assumption is the following: if the intrinsic dimension of the data is d, then the remaining d -D eigenvalues, corresponding to the eigenvectors that are discarded, should be significatively small. This principle is expressed using

Prop = d j=1 λ j / D j=1 λ j ,
which is equal to the proportion of the original variance kept. We will write W d the square matrix of size D containing the d eigenvectors corresponding of the higher eigenvalues, and all the other columns are null. Then thanks to the Eckart-Young theorem [START_REF] Wang | Geometric structure of high-dimensional data and dimensionality reduction[END_REF] (theorem 5.1) it is possible to quantify the error of reduction of dimension such as :

Err PCA = V -W T d V W d 2 F = D j=d+1 λ 2 j (1.8)

Multi Dimensional Scaling (MDS) -A linear method

Goal Multidimensional scaling [START_REF] Schiffman | Introduction to multidimensional scaling: Theory, methods, and applications[END_REF] is a data mining technique used to reduce the dimensionality of the data by keeping as close as possible the pairwise distance between the data point. More technically, we have:

F = {v i } n i=1 -→ F = {v i } n i=1 (1.9) with v i -v j v i -v j , ∀i, j ∈ [1, n] 2
, where v iv j represents the Euclidean distance between v i and v j . So the goal is to find a configuration of the dataset in a lower dimension such as the similarities are preserved. There are two main objective functions to use in the corresponding optimization problem. The simplest one is the raw stress function, given by:

Φ(F ) = i,j∈[1,n] 2 ( v i -v j 2 2 -v i -v j 2 2 ),
(1.10)

STATE-OF-THE-ART ON DIMENSIONALITY REDUCTION

The second objective function is the Sammon cost function :

Φ(F ) = 1 i,j∈[1,n] 2 v i -v j 2 i,j∈[1,n] 2 ( v i -v j 2 2 -v i -v j 2 2 ). (1.11)
The unique difference between these two cost functions comes from the term: i,j∈[1,n] 2 v iv j 2 that is a way to highlight data points on F that are more similar.

Mathematical rationale

In the following we still represent our initial data by : F = {v i } n i=1 ∈ R D , and we would like to reduce the dimensionality and to build the representation

F = {v i } n i=1 ∈ R d . First, let us introduce some basic notions. Definition 1 The Euclidean distance between two vectors a = [a 1 , . . . , a D ] ∈ R D and b = [b 1 , . . . , b D ] ∈ R D is just: d 2 (a, b) = a -b 2 = D i=1 (a i -b i ) 2 .
To simplify the writting of this section, we will omit the subscript 2 for the L 2 norm and write just • instead of • 2 . For the rest of this part we will denote the (symmetric) Euclidean distance matrix and the squared Euclidean distance matrix, respectively with matrix D and S such as:

D i,j = d 2 (v i , v j ), S i,j = d 2 (v i , v j ) 2 . Definition 2 The inner product between two vectors a = [a 1 , . . . , a D ] ∈ R D and b = [b 1 , . . . , b D ] ∈ R D is defined as: < a, b >= D i=1 a i × b i .
Since the basic idea of MDS is to transform the distance matrix into a cross product matrix, we are going to introduce the Gram matrix, also called cross product matrix, written G, such as:

G i,j =< v i , v j > .
Another definition of the Gram matrix is the following one.

Definition 3 Let X be a set of data, the Gram matrix G associated to X is defined by :

G = X T X Since d 2 (v i , v j ) = √ < v i -v j , v i -v j >, we have d 2 (v i , v j ) = < v i , v i > + < v j , v j > -2 < v i , v j > thus D i,j = G i,i + G j,j -2G i,j .
By performing a double centering, we can transform the squared distance matrix into a Gram matrix G. To do it, we need a "centralizing matrix" H such as:

H = I - 1 n • 1 1 T , (1.12)
with 1 ∈ M n,1 (R) and

1 T = [1, . . . , 1].
CHAPTER 1. STATE-OF-THE-ART ON STATISTICAL LEARNING APPLIED TO HYPERSPECTRAL IMAGES 19 Lemma 4 The centralizing matrix H satisfies the following properties:

(1) H 2 = H;

(2) 1 T H = H 1 = 0 ;

(3) If X is a centred dataset, which means that columnwise mean(X) = 0, then X H = X.

Proof.

(1)

H H = (I - 1 n • 1 1 T ) (I - 1 n • 1 1 T ), H H = I - 2 n • 1 1 T + 1 n 2 • (1 1 T ) × (1 1 T ), however 1 1 T 1 1 T = n • 1 1 T , so H H = I - 1 n • 1 1 T , H H = H.
(2)

1 T H = 1 T (I - 1 n • 1 1 T ), 1 T H = 1 T - 1 n • 1 T (1 1 T ), however 1 T 1 1 T = n • 1 T , so 1 T H = 0, since H T = H, finally H 1 = 0.
(3) If X is a centred dataset, we have X 1 1 T = 0. Then we obtain

X H = X - 1 n • X (1 1 T ).
Lemma 5 Let X be a centred dataset, and G its Gram matrix, then we have

H G H = G.
We add a superscript c on G and X to note that they are centred: G c and X c .

Proof. Thanks to the Lemma 4, one has X c H = X c , and using Definition 3, we can write

G c = (X c ) T X c = H (X c ) T X c × H = H G c H.
Theorem 6 Let X be a dataset, and D be its Euclidean distance matrix. We have

D i,j = G c i,i + G c j,j -2G c i,j . (1.13) 1.3. STATE-OF-THE-ART ON DIMENSIONALITY REDUCTION Proof. Let us consider X = {x i } n i=1 ∈ R k , X c = {x c i } n i=1 ∈ R k , D, G c
, receptively a set of data, the same set of data but centred, the Euclidean distance matrix of the set X, the centred Gram matrix of X, and k ∈ N

D i,j = d 2 (x i , x j ) = d 2 (x i -mean(X), x j -mean(X)) = d 2 (x c i , x c j ), then, one has: D i,j = < x c i -x c j , x c i -x c j >, D i,j = < x c i , x c i > -2• < x c i , x c j > + < x c j , x c j >, D i,j = G c i,i + G c j,j -2G c i,j .
Theorem 7 Let X be a dataset, G its Gram matrix, S its squared distance matrix, centred matrices G c and S c obey the following relationship:

G c = - 1 2 S c . (1.14)
Proof. Using Lemma 5, we have

H G c H = G c . Let us translate it matricially, using G c = {G c ij } i,j∈[1,n] 2 , to obtain: H G c H = I - 1 n • 1 1 T G c I - 1 n • 1 1 T , thus H G c H = G c - 1 n • 1 1 T G c - 1 n • G c 1 1 T + 1 n 2 • 1 1 T G c 1 1 T = G c , hence, 1 1 T G c + G c 1 1 T - 1 n • 1 1 T G c 1 1 T = 0.
On the other hand, we have:

1 1 T G c =    1 . . . 1 . . . . . . . . . 1 . . . 1    ×    G c 11 . . . G c 1n . . . . . . . . . G c n1 . . . G c nn    ,
Then by doing the product, we get:

1 1 T G c =    n i=1 G c i1 . . . n i=1 G c in . . . . . . . . . n i=1 G c i1 . . . n i=1 G c in    .
In the same way, we can write:

G c 1 1 T =    n j=1 G c 1j . . . n j=1 G c 1j . . . . . . . . . n j=1 G c nj . . . n j=1 G c nj    .
And also:

1 1 T G c × 1 1 T =    i,j∈[1,n] 2 G c ij . . . i,j∈[1,n] 2 G c ij . . . . . . . . . i,j∈[1,n] 2 G c ij . . . i,j∈[1,n] 2 G c ij    .
CHAPTER 1. STATE-OF-THE-ART ON STATISTICAL LEARNING APPLIED TO HYPERSPECTRAL IMAGES 21

So finally we have:

   n i=1 G c i1 . . . n i=1 G c in . . . . . . . . . n i=1 G c i1 . . . n i=1 G c in    +    n j=1 G c 1j . . . n j=1 G c 1j . . . . . . . . . n j=1 G c nj . . . n j=1 G c nj    - 1 n •    i,j∈[1,n] 2 G c ij . . . i,j∈[1,n] 2 G c ij . . . . . . . . . i,j∈[1,n] 2 G c ij . . . i,j∈[1,n] 2 G c ij    = 0 n (R).
By multiplying this expression on the right by 1 we obtain:

   i,j∈[1,n] 2 G c ij . . . i,j∈[1,n] 2 G c ij    + n •    n j=1 G c 1j . . . n j=1 G c nj    -    i,j∈[1,n] 2 G c ij . . . i,j∈[1,n] 2 G c ij    = 0 n,1 (R) Hence ∀i ∈ [1, n] we have n j=1 G c ij = 0 and thanks to the symmetry of G c , ∀j ∈ [1, n] we have n i=1 G c ij = 0. In addition, it is known now that D 2 i,j = G c i,i + G c j,j -2G c i,j
thanks to Theorem 6, so finally:

n i=1 D 2 i,j = n i=1 G c i,i + n • G c j,j , (1.15)
n j=1 D 2 i,j = n • G c i,i + n j=1 G c j,j . (1.16)
But we also have that

S c = H S H = I - 1 n • 1 1 T S I - 1 n • 1 1 T , S c = S - 1 n • 1 1 T S - 1 n • S 1 1 T - 1 n 2 • 1 1 T S 1 1 T .
By rewriting matricially, and using S = {D 2 ij } i,j∈[1,n] 2 and all the matrix operations that we have done just above, we have:

S c ij = D 2 ij - 1 n • n i=1 D 2 i,j - 1 n • n j=1 D 2 i,j + 1 n 2 • i,j∈[1,n] 2 D 2 ij .
(1.17)

Using (1.15) and putting those equations on (1.17), we obtain the following new expression:

S c ij = D 2 ij - 1 n •   n i=1 G c i,i + n • G c j,j -n • G c i,i + n j=1 G c j,j - 1 n • i,j∈[1,n] 2 D 2 ij   S c ij = D 2 ij -G c ii -G c jj
According to Theorem 6, we have as expected:

S c ij = -2 • G c ij . 1.3. STATE-OF-THE-ART ON DIMENSIONALITY REDUCTION Definition 8 The Frobenius norm of a matrix A = {a ij } i,j∈[1,n] 2 ∈ M n (R) is given by A F = i,j∈[1,n] 2 |a ij | 2 .
We remind a classical result related to the Frobenius norm.

Proposition 9 Let us consider A, U ∈ M n (R) 2 such as U F = 1, then we have:

A U F = U A F = A F .
Going back to the MDS, our goal is to find F that minimizes:

Φ(F ) = i,j∈[1,n] 2 ( v i -v j 2 -v i -v j 2 ).
Let us call F the matrix associated to the dataset

{v i } i∈[1,n] ∈ R D , F the matrix associated to {v i } i∈[1,n] ∈ R d , S
and S respectively the Euclidean squared distance matrices of F and F , and G and G respectively their Gram matrices. Another way to see the optimization problem of MDS is to seek for a matrix S that minimizes:

Φ(S ) = i,j∈[1,n] 2 (S ij -S ij ), with rank(S ) = d, (1.18)
such that

i,j∈[1,n] 2 (S ij -S ij ) ≤ i,j∈[1,n] 2 |S ij -S ij |.
Thanks to the norms equivalence relationship in M n (R), we have:

∃α ∈ R + , S -S 1 = i,j∈[1,n] 2 |S ij -S ij | ≤ α • S -S F
Therefore, our new goal is to find S that minimizes:

Φ(S ) = S -S F , with rank(S ) = d, (1.19)
Then, using the proposition 9 with H = ( H F ) -1 • H, we have:

S -S F = H (S -S) H F = 1 H 2 F • H (S -S) × H F , S -S F = 2 H 2 F • G -G F .
Finally, we conclude that our problem is equivalent to find G that minimizes:

Φ(G ) = G -G F , with rank(S ) = d. (1.20)
To do it, we use now the classical Eckart and Young theorem [START_REF] Eckart | The approximation of one matrix by another of lower rank[END_REF].

CHAPTER 1. STATE-OF-THE-ART ON STATISTICAL LEARNING APPLIED TO HYPERSPECTRAL IMAGES 23

Theorem 10 Let us consider a matrix A ∈ M nm (R), with (n, m) ∈ N 2 . We can write its singular value decomposition A = U A Σ A V T A , where Σ A =diag(σ 1 , . . . , σ r ) is diagonal and contains the singular value of A sorted in decreasing order, and U A ∈ M nr (R) and V A ∈ M rm (R) have orthogonal vectors that contain the left and right singular vector of A, and where r = rank(A). Then if we call A ∈ M nm (R) the best d-rank approximation of A under the Frobenius norm which solution to:

A -A F , with rank(S ) = d is A = U A,d Σ A,d V T A,d
, where Σ A,d is diagonal and contains the top d singular value of A sorted in decreasing order, and U A ∈ M nr (R) and V A ∈ M rm (R) contain respectively the associated left and right singular vector of A. Moreover the error of minimisation is given by:

A -A F = r i=d+1 σ 2 i .
Then, applying this theorem to our problem, we have

G = U G,d Σ G,d U T G,d , (1.21)
where G is squared, since G = (F ) T F , we finally have

F = Σ 1/2 G,d U T G,d , (1.22)
where Σ G,d is diagonal and contains the top d eigenvalue of G and U G ∈ M n (R) contains the associated eigenvectors.

Isomap -A nonlinear method

Let us consider again a set of n points F = {v i } n i=1 ∈ R D . Isometric mapping (Isomap) [START_REF] Tenenbaum | A global geometric framework for nonlinear dimensionality reduction[END_REF] is a generalization of MDS, where the pairwise distance between points (v i , v j ) is replaced by an estimation of the pairwise geodesic distance between points (v i , v j ). The basic idea consists in building a weighted graph which goal is to recover the local geometry of the manifold, each node of this graph is a point v i and each edge is an Euclidean distance between the nodes. Then, some edges are removed, in order to build local neighbourhoods. There are two ways to build the graph: either k-neighbourhood (i.e., the number of neighbours is fixed to k) or -neighbourhood (i.e., a threshold of value is applied on the distance). Dijkstra algorithm is used to compute the shortest paths in the graph. It turns out that under certain conditions regarding the sampling of the manifold, its convexity and the fact that the manifold is smooth, the graph distance happens to be a good approximation of the geodesic distance on the manifold. Therefore the cost function to minimize in Isomap is

Φ(F ) = i,j∈[1,n] 2 ( v i -v j 2 -d geo (v i , v j ) 2 ), (1.23)
where d geo is an approximation of the geodesic distance on the manifold.

STATE-OF-THE-ART ON DIMENSIONALITY REDUCTION

Locally linear Embedding (LLE) -A nonlinear method

Locally Linear Embedding (LLE) is an unsupervised learning algorithm introduced by Roweis and Saul [START_REF] Roweis | Nonlinear dimensionality reduction by locally linear embedding[END_REF] that reduces the dimensionality of a dataset of points F = {v i } n i=1 ⊂ R D based on their local structure. It is based on the intuition that each initial point ∀i ∈ [1, n] v i can be linearly reconstructed by its neighbours (k-neighbourhood or -neighbourhood). The error of reconstruction is measured given by:

Φ 1 (W ) = i v i - j W ij • v j 2 ,
where W ij represents the contribution of the jth coordinate over v i , with the constrain that j W ij = 1 and W ij = 0 if v j does not belong to the neighbourhood of v i . This constraint, together with the cost function, makes that the reconstructed data are invariant to rotations, rescaling, and translation. Then one transforms v i to v i thanks to a mapping consisting in translation, rotation, and rescaling. Since W ij are invariant to this mapping, they are also the weight of reconstruction of v i . Thus to find v i , we just need to minimize the following equivalent cost function:

Φ 2 (F ) = i v i - j W ij • v j 2 , (1.24)
with the constrain j F ij = 0, ∀j ∈ [1, n] and F T F = I.

Hessian Eigenmaps -A nonlinear method

Hessian Eigenmaps introduced by Donoho and Grimes [START_REF] Donoho | Hessian eigenmaps: Locally linear embedding techniques for high-dimensional data[END_REF] achieves an embedding by minimizing the Hessian functional estimated from the sampled manifold. One of the major advantages of Hessian Eigenmaps is that it can recover the true manifold geometry even if it is not convex. Let us consider a set of n points F = {v i } n i=1 ∈ R D that lie on a manifold F of intrinsic dimension d. For any point v i we have a local orthonormal coordinate system (v 1 i , . . . , v d i ) on the tangent space T v i (F). Let us consider the projection operator Pr Tv i (v j ), where v j is in a neighbourhood of v i . For any twice differentiable function g : F -→ R it is possible to define its tangent Hessian matrix as:

H tan g (v i ) k,l = ∂ ∂v k i ∂ ∂v l i g(Pr -1 Tv i (v i )).
Let us define the quadratic form

H(g) = F H tan g (v i ) 2
F dv which represents the average over the data manifold F of the Frobenius norm of the Hessian matrix of g. Donoho and Grimes [START_REF] Donoho | Hessian eigenmaps: Locally linear embedding techniques for high-dimensional data[END_REF] proposed and proved the following key result: H(g) has a d + 1 dimensional null-space which consists in the constant function and the d-dimensional function spanned by the isometric coordinates. From this result, it is possible to formulate an algorithm to approximate the computation of Hessian Eigenmaps.

Kernels and Kernel PCA -A nonlinear method

Kernel trick

Kernel methods are used to define nonlinear decision functions while using linear methods as basic ingredients. They allow to project the data into a nonlinear space CHAPTER 1. STATE-OF-THE-ART ON STATISTICAL LEARNING APPLIED TO HYPERSPECTRAL IMAGES 25 that can be even of infinite dimension. For more information on kernels the reader is referred to [START_REF] Schölkopf | Learning with kernels: support vector machines, regularization, optimization, and beyond[END_REF].

Definition 11 A kernel is a function k : X × X → C which is symmetric and hermitian.

However most of the time, we work with positive definite kernels.

Definition 12 k is called a positive definite kernel if ∀{x 1 , . . . , x n } ∈ X n and ∀{α 1 , . . . , α n } ∈ C n , the following non-negativity condition holds:

n i=1 n j=1 α i α * j k(x i , x j ) ≥ 0.
If we consider a dataset {x 1 , . . . , x n } ∈ X n , thanks to positive definite kernel, it is possible to define a similarity matrix:

[K] ij := k(x i , x j ),
which is positive semi-definite (PSD). Usually this matrix is also called the Gram matrix, or the kernel matrix.

Positive definite kernels are considered as a generalization of inner product on nonlinear spaces, and so can be used to replace the inner product in algorithms. In addition, the function k : X × X → R, defined by:

∀(i, j) ∈ [1, n] k(x i , x j ) =< x i , x j > R D , where < .., .. > R D is the inner product on R D , is just called the linear kernel.
Definition 13 A Hilbert space H is a vector space with a real or complex inner product space that is also a complete metric space with respect to the distance function induced by the inner product. That means that every Cauchy sequence in H has a limit in H.

We remind now the important Moore-Aronszajn theorem for kernel methods.

Theorem 14 k is a positive definite kernel if and only if there exists a Hilbert space H and a mapping φ : X → H such that k(x i , x j ) =< φ(x i ), φ(x j ) > H .

Definition 15 Let (H, < .., .. > H ) be a Hilbert space consisting of functions of X in C. The function k : X × X → C is the reproducing kernel of H, provided that the latter admits one if and only if:

• For any element x ∈ X , the function k(x, .) : t → k(x, t) belongs to H; • For all x ∈ X , g ∈ H, the reproducing property is verified: g(x) =< g, k(x, .) > H .
If a reproducing kernel exists, then H is called a reproducing kernel Hilbert space (RKHS).

It can be shown that every positive definite kernel is the reproducing kernel of at most one unique RKHS of functions from X to C. Reciprocally, if H is a RKHS, then it has a unique reproducing kernel. Moreover every reproducing kernel is a positive definite kernel. A proof of these theorems may be found in [START_REF] Rakotomamonjy | [END_REF]. From these properties, it results that we can replace φ(x i ) by k(x i , .) and therefore to obtain:

k(x i , x j ) =< φ(x i ), φ(x j ) > H =< k(x i , .)), k(x j , .) > H .
A final key result, named the Representer theorem.

STATE-OF-THE-ART ON DIMENSIONALITY REDUCTION

Theorem 16 Let X be a set endowed with a positive definite kernel k, and H k the corresponding RKHS, and x 1 , . . . , x n ⊂ X a finite set of points. Let Ψ : R n+1 → R be a function of n + 1 variables, strictly increasing with respect to the last variable. Then, any solution to the optimization problem:

min f ∈H k Ψ(g(x 1 ), . . . , g(x n ), g H k ), (1.25)
admits a representation of the form:

∀x ∈ X , g(x) = n i=1 α i k(x i , x), (1.26)
where

g H k = √ < g, g > H k .
The concept of "kernel trick" corresponds to the idea that for any algorithm that uses only the inner product between the input vectors, we can do this implicitly, in a Hilbert space called feature space φ, by replacing each inner product by an evaluation of the kernel. This evaluation can be done without having a representation of the feature map. The kernel trick allows to operate in a space larger without having to explicitly calculate the coordinates of the data in this space.

Kernel PCA

Let us consider a set of n vectors

v i ∈ R D , ∀i ∈ [1, n].
Let us map our data onto another space H, that have some interesting properties, by the mapping:

φ = R D → H v i → φ(v i ) (1.27)
where φ is a function that may be nonlinear, and depends on the kernel. The choice of the kernel is a difficult question that will be discussed latter. Now let us suppose that we have chosen a kernel. The goal of the kernel PCA (KPCA) is to find the set of vectors {w j , j ∈ [1, D]} that minimize the following cost function:

min 1 n n i φ(v i )-< φ(v i ), w j > H k •w j 2 H k , ∀j ∈ [1, P ], (1.28)
where w j , j ∈ [1, n] are the so-called principal components. By developing the distance one gets:

φ(x i )-< φ(v i ), w j > H k •w j 2 H k = φ(v i ) 2 H k -2 < φ(v i ), w j > 2 H k + < φ(v i ), w j > 2 H k w j 2 H k (1.29)
By adding the constraint that w j 2 H k = 1, replacing (1.29) in 1.28 and keeping the term that depend on w j , we have the following new objective function:

max( 1 n × n i < φ(v i ), w j > 2 H k ), with w j 2 H k = 1, ∀j ∈ [1, P ].
(1.30) CHAPTER 1. STATE-OF-THE-ART ON STATISTICAL LEARNING APPLIED TO HYPERSPECTRAL IMAGES 27 Problem (1.30) can be rewritten thanks to the Lagrange multiplier theorem as follows

L(w j , λ) = 1 n n i < φ(v i ), w j > 2 H k -λ ( w j 2 H k -1) (1.31)
where λ ∈ R. Thanks to the Representer theorem presented above, w j can be written as:

w j = n l=1 α l,j φ(v l ). (1.32)
Therefore, finally, we have:

L(α j , λ) = 1 n n i n l=1 α l,j < φ(v i ), φ(v l ) > H k 2 - λ n l=1 n t=1 α l,j α t,j k(v l , v t ) -1 ,
which can be rewriten in a matrix way as:

L(α j , λ) = 1 n α t j k 2 α j -λ • α t j k α j -1 . (1.33)
By deriving it, we obtain:

∂L ∂α j (α j , λ) = 2 n k 2 α j -2λ k α j = 0. (1.34)
In conclusion, KPCA optimizaton problem is equivalent to eigenproblem:

k α j = nλ • α j (1.35)
such that (1.34) and (1.35) are equivalent for non null eigenvalue, in the other cases, the solutions are not interesting for the maximization problem, and would not be taken into consideration. Finally α j , j ∈ [1, D], are the eigenvectors of k, and the normalization impose that the eigenvector are:

αj = 1 λ j α j .

Importance of the choice of the kernel

The advantage of the kernel trick is that we can use many different kernels without having to compute explicitly the mapping φ(v i ). Two of the most popular kernels are:

• The polynomial kernel:

k(v i , v j ) = (< v i , v j > R D +c) P ,
where P is the degree of the kernel and c a scalar constant;

1.3. STATE-OF-THE-ART ON DIMENSIONALITY REDUCTION

• The radial basis function (rbf) kernel, often called Gaussian kernel:

k(v i , v j ) = e -v i -v j 2 R D 2σ 2
, with scaling parameter σ. This kernel bring the data into a RKHS of infinite dimension. Kernel function k(v i , v j ) can be seen as an interaction function between v i and v j , that quantifies the similarity between two objects. Here we work on high dimensional spaces, where the notion of similarity is a key issue to address the curse of dimensionality. Usually we want a good similarity measure to be high for data that belong to the same cluster and small on the contrary. However since we do not have access to this information, this involves that the question of similarity in high dimensional spaces is far from being simple.

We have represented in Figures 1.5 and 1.6 the results of KPCA for different kernels. These results illustrate the importance of the choice of the kernel. By the way, the choice of σ in the case of the rbf kernel has a strong impact on the dimensionality reduction, as shown in the figures.

In this thesis, and when a particular approach is not mentioned, we have chosen for the rbf kernel parameter

σ = median({ v i -v j R D , (i, j) ∈ [1, n] 2 }).
For readers interested on the choice of σ are referred to [START_REF] Takeuchi | Nonparametric quantile estimation[END_REF].

Dictionary learning

The goal of dictionary learning is to find a sparse representation which approximates an image from a set of consistent images. Let us consider that we have a collection of D images, denoted here by y j , with j ∈ [1, D]. Our goal is to find a dictionary D, composed of d atoms φ k , with k ∈ [1, d], such as each image y j can be expressed as a linear combination of atoms from the dictionary:

y j = k α j,k φ k = α T j φ, ∀j. (1.36)
In order to find a unique dictionary, some constraints of sparsity are imposed on α j . Typically the objective is to find a sparse vector α j that would contain a small number of non-zero coefficients. The reformulated optimization problem is written as follows [START_REF] Tosic | Dictionary learning[END_REF]:

min    y j -α T j φ Term of estimation +λ α j 0 Term of sparseness    , ∀j ∈ [1, D].
(1.37)

In order to solve this problem, the proposed techniques can be divided into two families. First, supervised techniques, where the dictionary D is given, and so the goal is only to find α j . We do not consider this kind of approaches. Second, unsupervised techniques, where D and α j are both estimated, which can be subdivided into three main kinds: i) probabilistic methods, ii) clustering-based methods, and iii) learning a dictionary with specific structure.

Maximum likelihood method -A probabilistic method

The maximum likelihood method for dictionary learning [START_REF] Mairal | Supervised dictionary learning[END_REF][START_REF] Lewicki | Learning overcomplete representations[END_REF] leans on the hypothesis that the signals can be written as:

y j = α T j φ + n j , ∀j,
where n j is a white Gaussian residual vector with covariance σ 2 . Then the goal of the approach is to find α and D maximizing the likelihood P (Y |D) = Π D i=1 P (y i |D) with: P (y i |D) = P (y i , α|D)dα = P (y i |α, D)P (α)dα.

VQ objective -A clustering based method

In this technique [START_REF] Schmid-Saugeon | Dictionary design for matching pursuit and application to motion-compensated video coding[END_REF], where the datapoints are partitioned into K clusters with for instance K-means algorithm, each data cluster is approximated thanks to one atom of the dictionary, which usually corresponds to the central point (i.e., mean) of the cluster.

The K-SVD -A clustering based method K-SVD [START_REF] Aharon | An algorithm for designing overcomplete dictionaries for sparse representation[END_REF] is a powerfull generalization of VQ objective. It is contrary to most of the methods where the optimisation problem (1.37) is solved thanks to a double optimization problem, where one first fix α, and find D, and then, one fixes D, to find α. In K-SVD each atom φ k of D is update sequentially using an Singular Value Decomposition (SVD), then thanks to the SVD, a new atom φk is found as well as its corresonding coefficients. This technique is in a way a generalisation of the VQ objective since each data y j can be expressed thanks to multiple atoms, from multiple partitions.

Unions of orthonormal Bases -Structured dictionary

The technique [START_REF] Lesage | Learning unions of orthonormal bases with thresholded singular value decomposition[END_REF] is based on a particular constrained optimization problem, considering that the solution is a set of union of orthonormal bases (1.37). Its rationale is founded on the fact that many signals can be seen as a set of orthonormal bases, and also that, thanks to this structure, the problem of optimization is easier to solve. In this method, each bases is updated sequentially, but the optimization problem does not update the atom of the basis and the corresponding coefficients at the same time.

State-of-the-art on classification 1.4.1 Regression

Supervised classification can be seen as a special case of regression. For this reason, let us start by introducing a background on regression [START_REF] Bishop | Pattern recognition[END_REF].

Let us suppose that we have n datapoints (or vectors)

{v i } ∈ R D , with i ∈ [1, n].
Each of these points is associated to a target value {t i } ∈ R. The set {v i , t i }, i ∈ [1, n], is called the training set. The goal of regression is to predict the target value t i for a new data v i with i / ∈ [1, n]. To be able to predict such a value t i , we can consider that it can be written as a linear combination of the coordinates of v i . Let us write v i = (v i,1 , . . . , v i,j , . . . , v i,D ) t ∈ R D . A simple model often used in regression is to consider that the prediction function is given by:

y(ω, v i ) = ω 0 + ω 1 v i,1 + . . . + ω D v i,D = ω 0 + D j=1 ω j v i,j .
(1.38)

Our goal is to learn the parameters {ω 0 , . . . , ω D } thanks to the training set. This model (1.38) is called linear regression, and may have some limitations. That is why we prefer to consider a more general model:

y(ω, v i ) = ω 0 + M j=1 ω j φ j (v i ), (1.39)
where φ j (v i ) represents a set of basis functions, with M the number of basis. There are many choices of basis functions. For instance, we can use a polynomial basis

φ j (v i ) = v k i,j , (1.40)
where the power k is a hyper-parameter or we can consider the rbf:

φ j (v i ) = exp - v i -µ j 2 2σ 2 , (1.41)
where now the hyper-parameters are µ j , which governs the spatial locations, and σ the scaling factor. Let us consider that the target data is given by the previous deterministic function, corrupted by Gaussian noise of zero mean Gaussian and inverse variance β, such that:

t i = y(ω, v i ) + ,
where y(ω, v i ) = ω t φ(v i ), with ω t = (ω 0 , . . . , ω M ), φ(v i ) = (φ 0 (v i ), . . . , φ M (v i )) t and φ 0 (v i ) = 1. Let us denote E the random variable associated to the noise such that E ∼ N (0, β -1 ). We can write τ for the random variable associated to the target value, such that we have τ ∼ N (y(ω, v i ), β -1 ), which depends on two parameters, ω and β and a spatial position v i on the manifold of the data. One can see that this is related to the notion a random functions that will be introduce in the sequel of the thesis. Let us consider that the training set is drawn independently from the previous law. Then we can write the likelihood function of the parameters ω and β:

L(t 1 , . . . , t n /ω, β) = n i=1 √ β √ 2π exp -β(t i -y((ω, v i )) 2 2 .
Taking the logarithm of the likelihood function, we have:

log L (t 1 , . . . , t n /ω, β) = n i=1 1/2. log β -1/2 log 2π -β/2(t i -y((ω, v i )) 2 .

STATE-OF-THE-ART ON CLASSIFICATION

If we want to find the set of parameters that maximize the likelihood, we have first to derive it according to each of the parameters of the log-likelihood, and equals it to zero. On the previous expression the term that depends just on ω is:

E d (ω) = β 2 n i=1 (t i -y((ω, v i )) 2 .
Our goal is to find ω that minimizes E d . One can recognize the typical cost function of machine learning regression. In order to minimize this cost function, we can derive it, equal it to zero, to finally obtain that:

ω M L = (Φ t Φ) -1 Φ t t, (1.42)
where Φ ∈ M n,M +1 (R) is defined by

Φ =    φ 0 (v 1 ) . . . φ M (v 1 ) . . . . . . . . . φ 0 (v n ) . . . φ M (v n )    (1.43)
and where t = (t 1 , . . . , t n ) t is the vector of all the training targets values. It is also possible to estimate β M L as:

β M L = 1 n n i=1 t n -ω t M L φ(x i ) 2 , (1.44)
such that β M L provides us information on the precision of the regression. This is the classical case of regression. However if one has access to some prior knowledge on ω, this can help to regularize the cost function. Let us consider that the random variable associated with ω follows a normal distribution N (m 0 , S 0 ), which is called the prior distribution. Then we need to compute the posterior distribution and we can estimate the new value of ω. As shown in [START_REF] Bishop | Pattern recognition[END_REF] (page 153), in the case where m 0 = 0 ∈ R M +1 and S 0 = α • I, with I the identity matrix of size M + 1, then the cost function is called ridge regression is just given by:

E B (ω) = β 2n n i=1 (t i -y((ω, v i )) 2 + α • ω t ω, (1.45)
whose solution is equal to:

ω B = Φ t Φ + α β • I -1 Φ t t (1.46)
where this new term α/β • I can be seen as a way to improve the inversion of the matrix, since small eigenvalues that might be problematic are improved thanks to this term. In addition, in the regularization theory [START_REF] Friedman | The elements of statistical learning[END_REF][START_REF] Schölkopf | Learning with kernels: support vector machines, regularization, optimization, and beyond[END_REF], it can also be seen as a way to avoid overfitting on the training set.

In the previous case, we were looking for a ω ∈ R M +1 of finite dimension. By means of the kernel trick, it is now possible to have a potential infinite dimension ω. Let us map the training set v i ∈ R D , ∀i ∈ [1, n], onto another space using

φ = R D → H x i → φ(x i ) (1.47)
CHAPTER 1. STATE-OF-THE-ART ON STATISTICAL LEARNING APPLIED TO HYPERSPECTRAL IMAGES 33 where again, φ is a function that may be nonlinear, and depends on the choice of the kernel. Then, in our regression case, the cost function can be now reformulated as:

E B (ω) = 1 2n n i=1 (t i -y(ω, v i )) 2 + λ • ω t ω. (1.48)
Apparently nothing is different, however since H might be infinite dimensional, and since we might not have access to φ, then the problem has changed. In order to solve this problem, we use the Representer theorem, and say that a solution of the problem must be of the form:

y(ω, v p ) = n i=1 γ i k(v p , v i ), (1.49)
where k is the kernel associated to φ, and where v p is not necessarily in the training set. At first v p would be in the training set to learn the parameters γ i , then it will not be. We can rewrite this problem of regression in a matrix way:

E B (γ) = 1 2n (Kγ -t) t (Kγ -t) + λ • γ t Kγ, (1.50)
where γ = (γ 1 , . . . , γ n ) and K is the Gram matrix associated with the kernel such that K ij = k(v i , v j ). The expression (1.50) can be developped as follows:

E B (γ) = 1 2n γ t K 2 γ -2t t Kγ + t t t + λ • γ t Kγ.
Then, by deriving it and setting equal to zero we obtain:

γ = (K + λ • I) -1 t, (1.51)
and finally, the solution is:

y(ω, v p ) = γ t k(v p , :), (1.52)
with k(v p , :) = (k(v p , v 1 ), . . . , k(v p , v n )) t . It happens that in the case λ = 0, the solution of the kernel regression is the same as the one of the kriging interpolation without constrain in the coefficients of the linear combination, see corresponding chapter.

Classification From regression to classification

The problem of classification can be considered as a special case of regression, where the target value is assigned to one of the K discrete classes. Then the input space is divided into regions, where each region is assigned to one class. Thus the goal of classification algorithms is to find the decision boundaries between the classes. In classification, the target vector t i of for instance class 2 is typically written as:

t i = (0, 1, 0, . . . , 0) .
(1.53)

STATE-OF-THE-ART ON CLASSIFICATION

As discussed above, in linear regression, the model is of the form: y(v i ) = ω 0 + D j=1 ω j v i,j . It turns out that in classification, since most of the time we want that the result of y is in [0, 1], it is possible to use an "activation function" a such that the model becomes:

y(v i ) = a ω 0 + D j=1 ω j v i,j .
(1.54)

From binary hyperplane classification to multi-class hyperplane classification

Binary classification At first, we address the simple case of binary classification.

We denote these two classes t i ∈ {-1, +1}. In addition, we consider that

y(v i ) = ω 0 + D j=1 ω j v i,j = ω 0 + ω t v i .
In the rest of the chapter, we write b 0 = ω 0 , since it represents a bias. We look for the weights ω such that if the data point v i is of the class 1 then y(v i ) ≥ 0, and y(v i ) < 0 if the point is of class -1. Hence, this activation function is just sign(x), which gives the sign of x, thus it separates the space in two classes.

Let us consider two points v 1 and v 2 which are on the decision boundary, therefore y(v 1 ) = y(v 2 ) = 0. Thus ω t (v 1v 2 ) = 0, then ω is orthogonal to the decision boundary. Finding ω will therefore determine the direction of the decision boundary. Moreover if v 1 is a vector on the decision boundary, y(v 1 ) = 0, so

ω t v 1 ω = - b 0 ω .
This implies that the distance from the origin to the decision surface is determined by ω 0 . We illustrate this property in Figure 1.7, for a case of D = 2. In addition, for any point v i , we can decompose it into the base composed of the information tangent to the decision boundary v ⊥ i and the information tangent to the w, i.e.,

v i = v ⊥ i + r ω ω ,
where r ∈ R represents a coordinate. Then we have on the one hand:

ω t v i = ω t v ⊥ i + r ω ω ,
and on the other hand, we add ω 0 and use the fact that y(v ⊥ i ) = 0, to obtain that

y(v i ) = r ω ω . (1.55)
Instead of considering the classification in the linear space, one might work also on the kernel space. To do that, we just need a feature space φ and then to apply the classification on this space thanks to These issues can be avoided by considering a single K-class classification problem, and then to attribute to each target:

y(v i ) = ω 0 + ω t φ(v i ).
y(v i ) k = ω t k v i + b 0,k , (1.56)
such that each point is assigned to the class k such that y(v i ) k > y(v i ) j , for all j ∈ [1, K]. For more details about this technique see [START_REF] Bishop | Pattern recognition[END_REF]. More recent techniques have been invented that provide interesting results, see for instance [START_REF] Crammer | On the algorithmic implementation of multiclass kernel-based vector machines[END_REF]. We conclude this paragraph by saying that multiclass classification is a rather difficult; open problem; since there are many alternative algorithms to transform from binary to multiclass classification.

Support Vector Machine (SVM)

We have now all the ingredients required to understand the Support Vector Machine (SVM) classification procedure. We start with the same binary classification problem:

y(v i ) = b 0 + ω t φ(v i ).
Since the weights ω might be an infinite dimensional vector, we change previous notation to the kernel one: y(v i ) = g(v i ), with g being a function of the Hilbert space H associated with φ.

Since the distance of the data according to the hyperplane determines their classes, we might be interested in finding a hyperplane that is far away from the data. The corresponding distance is called the margin. According to the discussion of previous paragraphs, the distance of each point to the hyperplane is in binary classification given by

t i g(v i ) g H . (1.57)
Hence the solution to find the maximum marging to the closest point v i is the decision boundary, which can be found by solving:

arg max g∈H min i g(v i ) g H . (1.58)
The intuition behind this idea is that for the definition of the hyperplane we focus just on the vectors which are near the boundary of influence. This can be formulated thanks to the hinge loss function, i.e.,

Φ(t i , y(g, v i )) = max (0, 1 -t i g(v i )) . (1.59)
This function is equal to zero if the data are on the correct side of the hyperplane. For data on the wrong side of the decision boundary, the value of the function is proportional to the distance from the margin. Thus the empirical risk function to be minimized is:

min f ∈H 1 n n i=1 Φ(t i , y(g, v i )) + λ g H , (1.60)
where the parameter λ determines the trade-off between increasing the margin-size and ensuring that the v i lie on the correct side of the margin, and it help also to regularize the solution.

CHAPTER 1. STATE-OF-THE-ART ON STATISTICAL LEARNING APPLIED TO HYPERSPECTRAL IMAGES 37 Using one more time the Representer theorem, the solution of the previous problem can be expanded as:

g(v p ) = n i=1 α i k(v i , v p ).
(1.61)

By plugging this formulation into the original problem, the following convex optimization problem is obtained:

min α∈R n 1 n n i=1 Φ(t i n j=1 α i k(v i , v j )) + λ n i=1 n j=1 α i α j k(v i , v j ).
(1.62)

Different loss functions provide different classification algorithms. Here we use the hinge function. However, this function is not differentiable, so one should reformulate the problem as follows:

min α∈R n ,ξ∈R n 1 n n i=1 ξ i + λ n i=1 n j=1 α i α j k(v i , v j ), (1.63)
subject to:

t i n j=1 α i k(v i , v j ) + ξ i -1 ≥ 0 ξ i ≥ 0 (1.64)
The Lagrangian of this problem is given

L(α, ξ, µ, ν) = 1 n n i=1 ξ i + λ n i=1 n j=1 α i α j k(v i , v j ) + n i=1 µ i (t i n j=1 α i k(v i , v j ) + ξ i -1) + n i=1 ν i ξ i ,
which corresponds to a classical minimization of a convex quadratic function with linear constraints that can be solved using a quadratic program optimization package.

Neural network classification

The architecture (Artificial) neural networks are a classification approach which attempt to find a mathematical representation of how our biological system processes information. First, let us focus on the model. In classification, the optimization problem was modelled by:

y(v i ) = a(ω 0 + D j=1 ω j φ j (v i )), (1.65)
where φ is representing a feature space associated with a kernel k, where the kernel is fixed and not learned. Here the goal is to make the kernel depending on many parameters which provide a huge flexibility to the feature space. Moreover we would 1.4. STATE-OF-THE-ART ON CLASSIFICATION like that these parameters and, so the feature space, are learned during the classification process. Therefore, the classification will learn ω and φ. To construct the neural network, we first consider a linear model that depends on the input data:

c k = ω (1) 0,k + D j=1 ω (1) j,k v i,j .
(1.66)

with k ∈ [1, K] and where each c k is a neurone of the first layer. The superscript (1) indicates that these parameters are the parameters of the first layer. Then, a nonlinear activation function a is applied on these quantities c k :

z k = a (1) (c k ). (1.67)
We can choose different kinds of activation functions, typically:

• A sigmoid function: a = tanh; or

• Rectified Linear Unit (ReLU): a(x) = 0 if x < 0 x if x ≥ 0 .
On the output of the first layer, a second linear combination is applied:

d k1 = ω (2) 0,k1 + K k=1 ω (2)
k,k1 z k .

(1.68) with k1 ∈ [1, K1], which is followed by another activation function. We can combine the different stages and obtain:

y k1 (v i ) = a (2) ω (2) 0,k1 + K k=1 ω (2) k,k1 .a (1) ω (1) 0,k + D j=1 ω (1) j,k v i,j . (1.69)
This function can be represented in the form of a network diagram as shown in Figure 1.9. Then to evaluate the parameter ω we proceed to a forward and backward propagation of information through the network using typically the stochastic gradient descent, see the reference [START_REF] Bottou | Online learning and stochastic approximations[END_REF][START_REF] Lecun | Gradient-based learning applied to document recognition[END_REF].

The neural network model described above comprises just two stages of processing, where each step corresponds to a perceptron algorithm [START_REF] Bishop | Pattern recognition[END_REF][START_REF] Bishop | Neural networks for pattern recognition[END_REF]. Hence this neural network is referred as a multilayer perceptron (MLP). The first unit is called the input unit, the last is called the output unit, and between these two units we have what is called the hidden units. In the case where the activation functions of all the hidden units are just linear functions, the corresponding network is equivalent to a network without hidden units. This is due to the fact that the composition of linear operators gives a linear operator, so we just need to find one linear operator that will summarize all the hidden units. Moreover, if the number of hidden units is smaller than either the number of input or output units, the corresponding network simplifies the data by proceeding to a dimensionality reduction. Biological motivation As explained in [START_REF] Johnson | The stanford cs class cs231n: Convolutional neural networks for visual recognition[END_REF], the basic unit of the brain is a neuron, and our brain is composed of approximately 86 billion neurons. These neurons are connected thanks to the synapses, as illustrated in Figure 1.10. It has been inferred that each neuron receives an input signal from its dendrites and produces an output signal along one single axon. This axon can eventually be connected via synapses to the dendrites of different neurons. The basic idea is that the synaptic strengths play the role of the coefficient ω (l) i,j and they control the influence of the signal.

Figure 1.10: Illustration of the neural network from a biological point of view [START_REF] Johnson | The stanford cs class cs231n: Convolutional neural networks for visual recognition[END_REF].

In the basic models, the signal is carried to the dendrites; once it is received, all the signals are sent to the cell body where they all get summed up. If the final sum is superior to a certain threshold, the neuron is activated, sending a spike along its axon. This can be interpreted as an activation function.

Convolutional Neural Network Neural networks are difficult to optimize (to learn) since we consider all the coordinates of v i . That is why usually a sparse network that does not have all the possible connections between two layers is considered. One of this type of neural networks is underlying the paradigm of Convolutional Neural Networks (CNN). CNN are considered today as providing impressive improvement of state-of-the-are classification task in computer vision [START_REF] Krizhevsky | Imagenet classification with deep convolutional neural networks[END_REF]. Right 1.4. STATE-OF-THE-ART ON CLASSIFICATION now, one of the best results of classification are obtained with a CNN of 152 layers [START_REF] He | Deep residual learning for image recognition[END_REF]. In these networks, the layer is organized into planes which can be seen as a feature map. Each layer takes as input only data from a small spatial subregion of the image as represented in Figure 1.11. Let us use the notation of the "matconvnet library" [START_REF] Vedaldi | Matconvnet: Convolutional neural networks for matlab[END_REF] documentation in order to explain a bit more the CNN. An image is denoted as:

f : E -→ F (i, j) -→ f (i, j, :) ∈ R D (1.70)
such that if we want to access to a particular feature coordinate k at a spatial position (i, j) we write it f (i, j, k). A CNN can be seen as the composition of a given number of functions, each one of them corresponding to a layer:

y k1 (f ) = g N . . . g 1 f, ω (1) , ω (N ) (1.71)
Function g l , l = 1, • • • , N takes as input the data ot the previous layer and a set of parameters ω (l) which are learned from data in order to solve a target problem. There are different kinds of functions g l as illustrated in Figure 1.12. Some g l do not have any parameter. The first function is typically the convolution. We note that here it is a threedimensional convolution, in the sense that it operates on a spatial map taking into consideration all the feature channels:

y(i , j , k ) (1) = g f, ω (1) = i,j,k ω (1) i,j,k,k f (i + i -1, j + j -1, k).
(1.72)

with k ∈ [1, N 1 ],
and N 1 is the number of feature maps for layer 1. y(i , j , k

) (1)
represents the information of one feature maps k 1 of layer 1 in position (i , j ). ω (1) represents a multi-dimensional filter. As we said, a nonlinearity is required to obtain useful features. The simplest nonlinear layer is a ReLU:

y(i , j , k ) (2) = g y(i , j , k ) (1) = max(0, y(i , j , k ) (1)
).

(1.73)

Next typical layers are so-called pooling layer. A pooling operator works on individual feature channels, replacing nearby features values by one feature thanks to a suitable operator. The most popular are the max-pooling and the mean-pooling.

The max-pooling is defined by:

y(i, j, k) (3) = g y(:, :, k) (2) = max{y(i , j , k) (2) |i ≤ i ≤ i + p & j ≤ j ≤ j + p}, (1.74)
with p a scale parameter fixed by the user.

Chapter 2

Morphological Principal Component Analysis for Hyperspectral Image Analysis

Abstract This chapter deals with the issue of reducing the spectral dimension of a hyperspectral image using principal component analysis (PCA). To perform this dimensionality reduction, we propose to add spatial information in order to improve the features that are extracted. Several approaches proposed to add spatial information are discussed in this chapter. They are based on mathematical morphology operators. These morphological operators are the area opening/closing, granulometries and grey-scale distance function. We name the proposed family of techniques the Morphological Principal Component Analysis (MorphPCA). Present approaches provide new feature spaces able to handle jointly the spatial and spectral information of hyperspectral images. They are computationally simple since the key element is the computation of an empirical covariance matrix which integrates simultaneously both spatial and spectral information. The performance of the different feature spaces is assessed for different tasks in order to prove their practical interest.

Résumé

Ce chapitre traite de la réduction de la dimension spectrale d'une image hyperspectrale à l'aide de l'analyse des composantes principales (PCA). Pour réaliser cette réduction de dimensionnalité, nous proposons d'ajouter des informations spatiales. Plusieurs approches ont été proposées pour ajouter de l'information spatiale dans ce chapitre. Elles sont basées sur des opérateurs de morphologie mathématique.

Introduction

Hyperspectral images allow us to reconstruct the spectral profiles of objects imaged by the acquisition of several tens or hundred of narrow spectral bands. Conventionally, in many applications hyperspectral images are reduced in the spectral dimension before any processing. Most of hyperspectral image reduction methods are linear and do not care about the multiple sources of nonlinearity present in this kind 44 2.1. INTRODUCTION of images [START_REF] Guilfoyle | Further investigations into the use of linear and nonlinear mixing models for hyperspectral image analysis[END_REF]. Nonlinear reduction techniques are nowadays widely used on data reduction, and some of them have been used for hyperspectral images [START_REF] Bachmann | Exploiting manifold geometry in hyperspectral imagery[END_REF]. Nevertheless, most of these techniques present some disadvantages [START_REF] Van Der Maaten | Dimensionality reduction: A comparative review[END_REF] in comparison to the canonical linear principal component analysis (PCA). That is the rationale behind our choice of PCA as starting point. In particular, one major drawback of those nonlinear techniques is that they are computationally too complex in comparison to PCA. Hence most of the time, they cannot be applied on real full resolution images.

Another common disadvantage of both classical linear and nonlinear dimensionality reduction techniques is that they consider a hyperspectral image as a set of vectors.

They are appropriate when the data do not present useful spatial information, and therefore they are not totally adapted to images.

As mentioned above, dimensionality reduction in hyperspectral images is usually considered as a preprocessing step for supervised pixel classification as well as for other hyperspectral image tasks such as unmixing, target detection, etc. Hence, our goal is to incorporate spatial information into the dimensionality reduction (DR).

The contribution of our approach can be summarized as follows. We propose to add spatial information on the estimation of the covariance matrix used for PCA computation. This is done by means of morphological image representations, which involve a nonlinear embedding of the original hyperspectral image into a morphological feature space.

Many previous works have considered how to introduce spatial information into hyperspectral dimensionality image reduction. We can divide these techniques into different fields. The first family of techniques is close to our paradigm since they are based on mathematical morphology The rest of the chapter is organized as follows. Section 2 provides a remind on the mathematical morphology multi-scale representation tools used in our approach. Section 3 introduces in detail our approach named morphological principal component analysis (MorphPCA). In order to justify our framework, a summary of the classical theory underlying the standard PCA is provided as well as the notion of Pearson image correlation. Then, the four variants of MorphPCA are discussed, including an analysis of their corresponding covariance matrix meaning. The application of MorphPCA to hyperspectral dimensionality image reduction is considered in Section 4. That involves an assessment of the different variants according to different criteria. For some of the criteria, new techniques to evaluate the quality of dimensionality reduction techniques on image processing are introduced. Techniques arising from manifold learning are also considered in the comparaison. Finally, Section 5 closes the chapter with the conclusions.

The first data set was acquired over the city of Pavia (Italy), is a hyperspectral image of spatial size : 610 × 340 pixels, with 103 spectral bands. The second image, which represents the University of Houston, is a hyperspectral image of spatial size 349×1905 pixels and with 144 spectral bands [START_REF] Debes | Hyperspectral and lidar data fusion: Outcome of the 2013 grss data fusion contest[END_REF]. The last one called Indian Pines is a hyperspectral image of spatial size 145 × 145 pixels, and with 224 spectral bands. We note that this chapter is an extended and improved version of our following contributions [49,[START_REF] Franchi | Morphological principal component analysis for hyperspectral image analysis[END_REF]. CHAPTER 2. MORPHOLOGICAL PRINCIPAL COMPONENT ANALYSIS FOR HYPERSPECTRAL IMAGE ANALYSIS 45

Basics on morphological image representation

The goal of this section is to introduce a short background on morphological operators and transforms used in the sequel. Notation considered in the rest of this chapter is also stated.

Notation

Let E be a subset of the discrete space Z 2 , which represents the support space of a 2D image and F ⊆ R D be a set of pixel values in dimension D. Hence, it is assumed in our case that the value of a pixel x ∈ E is represented by a vector v ∈ F of dimension D, where discrete space E has a size of n 1 × n 2 pixels. This vector v represents the spectrum at position x. Additionally, we will write higher order tensors by calligraphic upper-case letters (I, S, . . .). The order of tensor

I ∈ R n 1 ×n 2 ×...×n J is J. Moreover if I ∈ R n 1 ×n 2 ×n 3 , for all i ∈ [1, n 3 ] I :,:,i represents a matrix of size n 1 × n 2
where the third component is equal to i. In our case we can also associate a tensor to the hyperspectral image F ∈ R n 1 ×n 2 ×D .

Nonlinear scale-spaces and morphological decomposition

Mathematical morphology is a well known nonlinear image processing methodology based on the application of complete lattice theory to spatial structures. Let f : E → Z be a grey-scale image. Area openings γ a s l (f ) (resp. area closings ϕ a s l (f )) are morphological filters that remove from the image f the bright (resp. dark) connected components having a surface area smaller than the parameter s l ∈ N [160]:

γ a s l (f ) = i {γ B i (f )|B i is connected and card(B i ) = s l }, (2.1) ϕ a s l (f ) = i {ϕ B i (f )|B i is connected and card(B i ) = s l }, (2.2)
where γ B (f ) and ϕ B (f ) represent respectively the morphological flat opening and closing according to structuring element B [START_REF] Serra | Image Analysis and Mathematical Morphology[END_REF]. We note that these connected filters can be implemented as binary filters on the stack decomposition of f into upper level sets. Figure 2.1 illustrates how area opening and area closing modify a simple image f . The image f in this toy example is composed of one black triangle of area equal to 30, 2 diamonds, one black and one white of area equal to 15. Finally the last connected components are 4 white circles and 5 black ones of area equal to 5. When an area opening is used (respectively closing) of threshold s l = 7, just the white (respectively black) circles are removed. Area opening and area closing are very relevant to simplify images, without deforming the contours of the objects remaining. In addition, area opening and closing can be used to produce a multi-scale decomposition of an image. The notion of morphological decomposition is related to the granulometry axiomatic [START_REF] Serra | Image Analysis and Mathematical Morphology[END_REF]. Let us consider {γ a s l }, 1 ≤ l ≤ S and {ϕ a s l }, 1 ≤ l ≤ S, two indexed families of area openings and closings respectively. Typically, the index l is associated to scale, or Figure 2.1: Illustration of an area opening γ a s l and an area closing ϕ a s l of image f , with s l = 7 pixels. We can see that the connected components removed by the opening operator are the white circles since their area is 5, so below 7, and similarly for the black circles in the closed image. more precisely to the surface area. Namely, we have on the one hand:

f = S l=1 (γ a s l-1 (f ) -γ a s l (f )) + γ a s S (f ); (2.3) f = ϕ a s S (f ) - S l=1 (ϕ a s l (f ) -ϕ a s l-1 (f )). (2.4)
On the other hand, we can rewrite the decomposition [START_REF] Velasco-Forero | Classification of hyperspectral images by tensor modeling and additive morphological decomposition[END_REF]:

f = 1/2 (γ s S (f ) + ϕ s S (f )) + S l=1 (γ a s l-1 (f ) -γ a s l (f )) - S l=1 (ϕ a s l (f ) -ϕ a s l-1 (f )) .
Therefore we have an additive decomposition of the initial image f into S scales, together with the average largest area opening and closing. We remark that the residue (γ a s l-1 (f )γ a s l (f )) represents bright details between levels s l and s l-1 . Similarly, (ϕ a s l (f )ϕ a s l-1 (f )) stands for dark details between levels s l and s l-1 . At this point, some issues must be taken into account. First, after decomposing an image into S scales, we have now to deal with an image representation of higher dimensionality. Second, the decomposition may not be optimal since it depends on the discretization of S scales, i.e., size of each scale. In order to illustrate that issue, we have represented in Figure 2.2(a) a channel of Pavia hyperspectral image and in Figure 2.2(b) its morphological decomposition by area openings that we have over-estimated. As it may be noticed from Figure 2.2(b), the choice of the scales is fundamental in order to avoid a redundant decomposition.

In order to deal with the problem of scale discretization, we propose to use the pattern spectrum that provides information about the image component size distribution. We can also notice another technique to find the optimal discretization [START_REF] Cavallaro | Automatic threshold selection for profiles of attribute filters based on granulometric characteristic functions[END_REF].

Pattern Spectrum

The notion of pattern spectrum (PS) [START_REF] Maragos | Pattern spectrum and multiscale shape representation. Pattern Analysis and Machine Intelligence[END_REF] corresponds to the probability density function (pdf) underlying a granulometric decomposition by morphological openings CHAPTER 2. MORPHOLOGICAL PRINCIPAL COMPONENT ANALYSIS FOR HYPERSPECTRAL IMAGE ANALYSIS 47 and closings [START_REF] Matheron | Random sets and integral geometry[END_REF][START_REF] Serra | Image Analysis and Mathematical Morphology[END_REF]. The area-based PS of f at size s l is given by

P S a (f, l) = Mes(γ a s l (f ) -γ a s l+1 (f )) /Mes(f ), (2.5) P S a (f, -l) = Mes(ϕ a s l+1 (f ) -ϕ a s l (f )) /Mes(f ), (2.6)
where Mes represents here the integral of the grey-scale image. Two images having the same pattern spectrum have the same morphological distribution according to the choice of the family of openings/closings. Since our goal is to have a nonredundant multi-scale representation with the same morphological representation than the original image, then by sampling the PS and choosing the scales of the distribution which keep it as similar as possible to the image PS, we can expect to find the appropriate discretization of scales. However, one can see in Figure 2.3, the PS is not a smooth function, and consequently, sampling it with a limited number of scales, would not lead to a good result.

Based on the analogy between the PS and probability density function, we can compute its corresponding cumulative pattern spectrum (CPS) for both sides l ≥ 0 and l ≤ 0. Naturally, this function is smoother than the PS. In order to select the appropriate scales, the CPS for openings and closings are sampled, where the number of samples is fixed and is equal to S, under the constraint that the sampled function must be as similar as possible to the original function.

An example of such sampling is given in Figure 2.3, where the approximation of the CPS is depicted in red and the CPS of the original image in blue. It is well known in probability that two distributions that have the same cumulative distribution function have the same probability distribution function. Based on this property, we can expect that the discretization from the CPS approximates the original PS of the image and consequently, the selected scales represent properly the size distribution of the image.

Grey-scale distance function

Let X be the closed set associated to a binary image. The distance function corresponding to set X gives at each point x ∈ X a positive number that depends on the position of x with respect to X and is given by [START_REF] Soille | Morphological image analysis: principles and applications[END_REF]:

dist(X)(x) = min{d(x, y) : y ∈ X c }, (2.7)
where d(x, y) is the Euclidean distance between points x and y, and where X c is the complement of set X. This well known transform is very useful in image processing [START_REF] Soille | Morphological image analysis: principles and applications[END_REF].

Distance function of binary images can be extended to grey-scale images f by considering its representation into upper level sets {X h (f )} a≤h≤b , where

X h (f ) = {x ∈ E : f (x) ≥ h} ,
such that a = min{f (x), x ∈ E}, and b = max{f (x), x ∈ E}. Then, the so-called grey-scale distance transform of f is defined as [START_REF] Molchanov | Distance transforms for real-valued functions[END_REF]:

dist(f )(x) = (b -a) -1 b h=a dist (X h (f )) (x).
(2.8)

That is, the grey-scale distance transform of f is equal to the sum of the distance functions from its upper level sets.

Morphological Principal Component Analysis

We introduce in this section the notion of Morphological Principal Component Analysis (MorphPCA) and its variants. Before reading this part we advice reader to read the section on the PCA the state of art chapter. On the next subsection we will focus on the covariance matrix needed to perform PCA.

Covariance matrix and Pearson correlation matrix

The covariance between two channels (or spectral bands) of an hyperspectral image F is computed as

Covar (F :,:,k , F :,:,k ) = 1 n n 1 i=1 n 2 j=1 [F i,j,k -E(F :,:,k )] [F i,j,k -E(F :,:,k )] , (2.9)
where E(F :,:,k ) is the mean of the hyperspectral image channel k. The covariance is very meaningful, however this is not a similarity measure [START_REF] Goshtasby | Similarity and dissimilarity measures[END_REF], in the sense of a metric, since it is not range limited. In order to fulfill this requirement, a solution consists in normalizing the covariance, which leads to the notion of Pearson correlation:

Corr (F :,:,k , F :,:,k ) = n 1 i=1 n 2 j=1 F i,j,k -E(F :,:,k ) σ k F i,j,k -E(F :,:,k ) σ k , (2.10)
where

σ k = 1 n n 1 i=1 n 2 j=1 (F i,j,k -E(F :,:,k )) 2 1/2
. The correlation coefficient varies between +1 and -1, such that Corr (F :,:,k , F :,:,k ) = 1 involves that F :,:,k and F :,:,k perfectly coincide. It has been proved that the best fitting case corresponds to [START_REF] Huntington | Mathematics and statistics, with an elementary account of the correlation coefficient and the correlation ratio[END_REF]:

F i,j,k = Corr (F :,:,k , F :,:,k ) σ k σ k (F i,j,k -E(F :,:,k )) + E(F :,:,k ). (2.11)
Therefore, from (2.11), we can see that the correlation is a linear coefficient between F i,j,k and F i,j,k . This means that Pearson correlation is a similarity criterion which depends on the intensities of the images and their linear relations.

MorphPCA and its variants

The fundamental idea of MorphPCA consists in replacing the covariance matrix V of PCA, which represents the statistical interaction of spectral bands, by a covariance matrix V Morpho computed from a morphological representation of the bands. Therefore, mathematical morphology is fully integrated in the dimensionality reduction problem by standard SVD computation to solve

V Morpho w j = λ j w j .
The corresponding principal components w j provide the projection space for the hyperspectral image F. This principle is illustrated in the diagram of

Scale-space Decomposition MorphPCA

In the first variant, we just use the area-based nonlinear scale-space discussed in the previous section. So the grey-scale image of each spectral band F :,:,k is decomposed into residues of area openings and area closings according to the discretization into S scales for each operator, i.e., r l (F :,:,k ) = γ a s l-1 (F :,:,k )γ a s l (F :,:,k ) and r -l (F :,:,k ) = ϕ a s l (F :,:,k )-ϕ a s l-1 (F :,:,k ), 1 ≤ l ≤ S. Thus we have increased the dimensionality of the initial dataset from a tensor (n 1 , n 2 , D) to a tensor (n 1 , n 2 , D, 2S + 1). As discussed in [START_REF] Velasco-Forero | Classification of hyperspectral images by tensor modeling and additive morphological decomposition[END_REF], this tensor can be reduced using high order-SVD techniques. We propose here to simply compute a covariance matrix as the sum of the covariance matrices from the various scales. More precisely, we introduce V Morpho-1 ∈ M D,D (R) with :

V Morpho-1 = S l=1 (V (l)) + S l=1 (V (-l))
(2.12)

where the covariance matrices at each scale l is obtained as

V (l) k,k = Covar (r l (F :,:,k ), r l (F :,:,k )) , 1 ≤ k, k ≤ D.
We note that it involves an assumption of independence of the various scales. We remark also that this technique is different from the classical approaches of differential profiles as [START_REF] Fauvel | Spectral and spatial classification of hyperspectral data using svms and morphological profiles[END_REF] where the morphological decomposition is applied after computing the spectral PCA (i.e., morphology plays a role for spatial/spectral classification but not for spatial/spectral dimensionality reduction as in our case).

Pattern Spectrum MorphPCA

In the second variant, we can consider a very compact representation of the morphological information associated to the area-based nonlinear scale-space of each spectral band. It simply involves considering the area-based PS of each spectral band as the variable to be used to find statistical redundancy on the data. In other words, the corresponding covariance matrix V Morpho-2 ∈ M D,D (R) is defined as :

V Morpho-2 k,k = Covar (P S a (F :,:,k , :), P S a (F :,:,k , :)) ,

(2.13) with 1 ≤ k, k ≤ D and where P S a (F :,:,k , l), -S ≤ l ≤ S, is the area-based pattern spectrum obtained by area-openings and area-closings. We note that the pattern spectrum can be seen as a kind of pdf of image structures. Consequently the Mor-phPCA associated to it explores the intrinsic dimensionality of sets of distributions instead of sets of vectors. For illustrating the information carried out by the PS, we have provided in Figure 2.9 the pattern spectra computed from three different bands of a hyperspectral image.

In order to better understand the interest of V Morpho-2 , we propose an analysis based on its Pearson correlation counterpart. Once the correlation of PS distribution is calculated, we have a linear coefficient between P S a (F :,:,k , l) and P S a (F :,:,k , l). However since the PS is the result of nonlinear operations, the underlying extracted features are naturally nonlinear.

Let us consider the two binary images of Figure 2.7(a), which represent two objects having exactly the same size. If the correlations are calculated, we have:

Corr (image1, image2) = 0 Corr (P S a (image1), P S a (image2)) = 1.
Hence, we can see that the morphological distribution being the same, the PS correlation is maximal. In a certain way, we observe that this transform builds sizeinvariants from the images and consequently it is robust to some groups of transforms and deformations. For instance, it is invariant to rotation and to translation.

Classical PCA on the spectral bands and the MorphPCA based on the PS can be compared by the corresponding correlation matrices from a hyperspectral image, such as the example plotted in Figure 2.10(a) and (b). From this visualization, we already observe that the bands are better discriminated between them.

Distance Function MorphPCA

Classical PCA for hyperspectral images is based on exploring covariances between spectral intensities. The previous MorphPCA involves changing the covariance into a morphological scale-space representation of the images. An alternative is founded when transforming each spectral band from an intensity based map to a metric based map where at each pixel the value is associated to both the initial intensity and the spatial relationships between the image structures. This objective can be achieved using the Molchanov grey-scale distance function [START_REF] Molchanov | Distance transforms for real-valued functions[END_REF] for each spectral band dist(F :,:,k ). The new covariance matrix V Morpho-3 ∈ M D,D (R) is now defined as:

V Morpho-3 k,k = Covar (dist(F :,:,k ), dist(F :,:,k )) , (2.14)
with 1 ≤ k, k ≤ D. Figure 2.9 depicts the corresponding grey-scale distance function from three spectral band of a hyperspectral image. We note that this function carries out simultaneously both intensity and shape information from the image. Let consider in detail the expression of the covariance of distance functions:

Covar (dist(F :,:,k ), dist(F :,:,k )) = Covar   max(F :,:,k ) h=min(F :,:,k ) d(X s (F :,:,k )), max(F :,:,k ) h =min(F :,:,k ) d(X h (F :,:,k ))   = max(F :,:,k ) h=min(F :,:,k ) max(F :,:,k ) h =min(F :,:,k ) Covar (d(X h (F :,:,k )), d(X h (F :,:,k ))) ,
where X h (F :,:,k ) denotes an upper level set at threshold h. The central term is the CHAPTER 2. MORPHOLOGICAL PRINCIPAL COMPONENT ANALYSIS FOR HYPERSPECTRAL IMAGE ANALYSIS 55 covariance between two binary distance functions and can be developed as follows:

Covar (d(X h (F :,:,k )), d(X h (F :,:,k ))) = E (d(X h (F :,:,k )), d(X h (F :,:,k ))) -E (d(X h (F :,:,k ))) E (d(X h (F :,:,k ))) = n -1 < d(X h (F :,:,k )), d(X h (F :,:,k )) > L 2 -E (d(X h (F :,:,k ))) E (d(X h (F :,:,k ))) ,
where < •, • > L 2 denotes the L 2 inner product. Using the classical relationship:

A -B 2 L 2 = A 2 L 2 + B 2 L 2 -2 < A, B > L 2 , ∀A, B ∈ R n ,
we finally obtain that:

Covar (d(X h (F :,:,k )), d(X h (F :,:,k ))) = (2n) -1 ( d(X h (F :,:,k )) 2 L 2 + d(X h (F :,:,k )) 2 L 2 ) -(2n) -1 d(X h (F :,:,k )) -d(X h (F :,:,k )) 2 L 2 -n -1 d(X h (F :,:,k )) 2 L 2 d(X h (F :,:,k )) 2 L 2 .
From this latter expression, the term

d(X h (F :,:,k )) -d(X h (F :,:,k )) 2 L 2 ,
can be identified as the Baddeley distance [START_REF] Baddeley | Errors in binary images and an lp version of the hausdorff metric[END_REF] used in shape analysis. This distance is somehow equivalent to the most classical Hausdorff distance between the upper level sets h of spectral band k and h of spectral band k . Thus, the underlying similarity from this covariance compares the shape of the spectral channels, and extracts a richer description than Pearson correlation from the spectral channels themselves. We note that the use of Hausdorff distance between upper level sets of hyperspectral bands was previously used in [START_REF] Velasco-Forero | Morphological image distances for hyperspectral dimensionality exploration using kernel-pca and isomap[END_REF]. Finally, to illustrate qualitatively the behavior of the distance function correlation, let us consider this time the three binary images depicted in Figure 2.7(b), where image 2 and image 3 represent the same object placed at a different location on the image. One has:

Corr (image1, image2) = Corr (image1, image3) Corr (dist(image1), dist(image2)) = Corr (dist(image1), dist(image3)) .
That is, this similarity criterion related to the use of distance function is more discriminative to the relative position of the objects on the image than the classical Pearson Correlation.

From Figure 2.10(c), one can compare now the correlation matrix using the greyscale distance function with the usual correlation matrix Figure 2.10(a). We note that this matrix provided also a better discrimination of bands cluster than the Pearson correlation matrix used in standard PCA.

Spatial/Spectral MorphPCA

As we have discussed, V Morpho-2 represents a compact morphological representation of the image, however the spectral intensity information is also important for dimensionality reduction. To come with a last variant of MorphPCA, we build another covariance matrix V Morpho-4 that represents the spectral and spatial information without increasing the dimensionality by the sum of two covariance matrices:

V Morpho-4 β = (1 -β)V + βV Morpho-2 ,
(2.15)

with β ∈ [0, 1], and where obviously V k,k = Covar (F :,:,k , F :,:,k ) and β stands for a regularization term that balances the spatial over the spectral information. This kind of linear combination of covariance matrices is similar to the one used in the combination of kernels, where kernels providing different information sources are combined to have a new kernel which integrates the various contributions [START_REF] Camps-Valls | Composite kernels for hyperspectral image classification[END_REF].

CHAPTER 2. MORPHOLOGICAL PRINCIPAL COMPONENT ANALYSIS FOR HYPERSPECTRAL IMAGE ANALYSIS 57 We can now use PCA and the four variants of MorphPCA to achieve dimensionality reduction (DR) of hyperspectral images. In order to evaluate the interest for such a purpose, it is necessary to establish quantitative criteria that should be assessed. These criteria will evaluate both locally and globally the effectiveness of the dimension reduction techniques.

(d) (e) (f) (g) (h) (i)
• Local criteria.

Criterion 1 (C1) The reconstructed hyperspectral image F using the first d principal components should be a regularized version of F in order to be more spatially sparse.

Criterion 2 (C2) The reconstructed hyperspectral image F using the first d principal components should preserve local homogeneity and be coherent with the original hyperspectral image F.

Criterion 3 (C3) The manifold of variables (i.e., intrinsic geometry) from the reconstructed hyperspectral image F should be as similar as possible to the manifold from original hyperspectral image F.

• Global criteria.

Criterion 4 (C4)

The number of bands d of the reduced hyperspectral image should be reduced as much as possible. It means that a spectrally sparse image is obtained.

Criterion 5 (C5) The reconstructed hyperspectral image F using the first d principal components should preserve the global similarity with the original hyperspectral image F. Or in other words, it should be a good noise-free approximation.

Criterion 6 (C6) Separability of spectral classes should be improved in the dimensionality reduced space. That involves in particular a better pixel classification.

These criteria are used to analyse the effectiveness of the DR methods studying locally and globally their ability to remove redundancy and to preserve the fully richness of the spectral and spatial information.

In order to assess C1, we compute the watershed transform [START_REF] Soille | Morphological image analysis: principles and applications[END_REF] on each channel F k of the hyperspectral image. Watershed transform is a morphological image segmentation approach which in a certain way can be seen as an unsupervised classification technique. The advantage of using the watershed is that it allows us to cluster the image according to the local homogeneity; thus, an image with less details will have less spatial classes than an image with many insignificant details. Then, the number of clusters N k of the watershed transform of F k is considered as an estimation of the image complexity. To evaluate the complexity of the reconstructed hyperspectral image, the number of spatial classes is counted after having done a watershed on each band. Finally, the mean of the number of spatial classes is taken, i.e.,

Error sparse spatially = (D -1 )

D k=1 N k .
Assessment of C2, which involves image homogeneity, is based on a partition of the image into homogenous regions. Let us first remind the definition of a α-flat zones [START_REF] Meyer | The levelings[END_REF], used for such a purpose. Given a distance d : R D ×R D -→ R, two pixels (f (x), f (y)) ∈ (R D ) 2 , from a vector-valued image f , belong to the same α-flat zone of f if and only if there is a path (p 0 , . . . , p n ) ∈ E n such as p 0 = x and p n = y and ∀i ∈ [1, n -1], d (f (p i ), f (p i+1 )) ≤ α, with α ∈ R + . Computing the α-flat zones for a given value of α produces therefore a spatial partition of the image into classes such that in each connected class the image values are linked by paths of local bounded variation. Working on the d eigenvectors, the image partition π α associated to the α-flat zones quantize spatially and spectrally an hyperspectral image, see example given in Figure 2.11. The goal of simultaneous spatial and spectral quantization of a hyperspectral image has been studied in [START_REF] Franchi | Morphological principal component analysis for hyperspectral image analysis[END_REF], where we have studied in detail the dependency on the distance. Moreover, we have shown that in high dimensional spaces quantization results are generally not good. For the case considered here, we propose to use the Euclidean distance on the reduced space by PCA or MorphPCA. The choice of α is done in order to guarantee a number C of α-flat zones similar for all the compared approaches. We can expect that, by fixing the number of zones in the partition, the difference between a partition and another one depends exclusively on the homogeneity of the image. Now, using the partition π α , the spectral mean value of pixels from the original image F in each spatial zone is computed. This quantization produces a simplified hyperspectral image, denoted F πα . Finally, we assess how far pixels of the original image from each α-flat zone are from their mean; which involves computing the following error

Error Homg = D k=1 n 1 ,n 2 i,j=1 |F i,j,k -F πα i,j,k | 2 .
This criterion can consequently be seen as a way to see the trustworthiness of the CHAPTER 2. MORPHOLOGICAL PRINCIPAL COMPONENT ANALYSIS FOR HYPERSPECTRAL IMAGE ANALYSIS 61

DR technique, since it measures if the homogeneous partition of the reduced hyperspectral image corresponds to the homogeneous zone of the original image. C3 has been evaluated by means of two manifold learning criteria called the Kintrusion and K-extrusion [START_REF] Lee | Quality assessment of dimensionality reduction: Rank-based criteria[END_REF]. They are based on other criteria called continuity and trustworthiness [START_REF] Venna | Local multidimensional scaling[END_REF]. These criteria reveal DR behavior in terms of its ability to preserve the data manifold structure. We have first sampled randomly 10 thousands spectra from our hyperspectral images, where each spectrum is a vector of dimension D. Then we have modelled the manifold by a graph where each node is a vector and each edge is the pairwise distance. We used the Euclidean distance as the pairwise distance. For the rest of the paragraph we note by x i a point from the original manifold, ν K i its neighbourhood of size K, xi the same point from the manifold after a DR and νK i its corresponding neighbourhood of size K. A neighbourhood of size K at point x i is composed of the K closest points to x i according to used metric. More precisely, the goal of K-extrusion is to measure how the points that were in the K-neighbourhood of x i are not preserved in the K-neighbourhood of xi after DR. The K-intrusion evaluates if points on the K-neighbourhood of xi on the DR manifold were in the K-neighbourhood of x i , i.e.,

M intrusion (K) = 1 - 2 G(K) × n i=1 j∈ν K i \ν K i r(i, j) -K, (2.16) M extrusion (K) = 1 - 2 G(K) × n i=1 j∈ν K i \ν K i r(i, j) -K, (2.17)
where r(i, j) is the rank of the data x j in the ordering according to the distance from x i , and respectively r(i, j) the rank of xj in the ordering according to the distance from xi , and the term G K scales the measure to be between zero and one, i.e.,

G(K) = N K(2N -3K -1) if K < N \ 2 N (N -K)(N -K -1) if K ≥ N \ 2.
(2.18)

For a better understanding of these formulae, see [START_REF] Venna | Neighborhood preservation in nonlinear projection methods: An experimental study[END_REF]. An important point is the dependence of these parameters on the size of the neighbourhood. From (2. [START_REF] Bishop | Pattern recognition[END_REF]) and (2.17), the following parameters are computed [START_REF] Venna | Neighborhood preservation in nonlinear projection methods: An experimental study[END_REF]:

Q(K) = M extrusion (K) + M intrusion (K) 2 , (2.19) B(K) = M intrusion (K) -M extrusion (K). (2.20)
The interest of Q(K) is that it estimates in average the quality of a DR technique, whereas B(K) reveals its behavior as being more intrusive or extrusive.

In order to assess C4, as classically done, the fraction of explained covariance is fixed. Then, the number of principal components needed is counted. The rationale is based on the fact that a good DR technique should reduce the number of dimensions and extract a limited number of features that would explain most of the image. However since this criterion is linked to a sparsity criterion, we would like to add a distortion criterion, C5.

The evaluation of C5 is founded on computing a pattern spectrum of both the original hyperspectral image and the DR image. An important point is that the 2.4. MORPHPCA APPLIED TO HYPERSPECTRAL IMAGES pattern spectrum will be computed by openings on the hyperspectral image viewed as a 3D image. By doing such assumption, the 3D openings are decomposing in a simultaneous way the spatial/spectral object of the image and the corresponding curves of the PS will represent the distribution of both the spatial and the spectral objects. Two hyperspectral images are similar if they have the same spectral/spatial size distribution. As discussed in Section 2, we prefer to use the cumulative PS in order to obtain a smoother curve. Normally we cannot deal with both spatial/spectral distortions with the reconstruction error of the two images. However we will also assess the SNR of the reconstruction error as an additional parameter.

Finally, C6 is related to supervised pixel classification of the hyperspectral image. We have considered the least square SVM algorithm [START_REF] Camps-Valls | Kernel-based methods for hyperspectral image classification[END_REF] as a learning technique, with a linear kernel or rbf kernel, where the rbf kernel is initialized for each DR technique using cross validation. For each supervised classification run, we used for the AVIRIS Indian Pine Image 5% of the available data as a training set and the remaining 95% to validate. For the ROSIS Pavia University image we use a subset of 50 spectra (about 1% of the available data) per class as a training set and the remaining spectra to validate.

Evaluation of algorithms

The studied DR techniques presented are listed and compared upon three mathematical and computing properties in Table 2.1. These properties were also considered in the excellent comparative review [START_REF] Van Der Maaten | Dimensionality reduction: A comparative review[END_REF]. For comparison, we have also included in the table the Kernel-PCA (KPCA), which is a powerful generalization of PCA allowing integrating morphological and spatial features into DR.

The first one is the number of free parameters to be chosen. The interest of having these free parameters is that it provides more flexibility to the techniques, whereas the related inconvenient is the difficulty for properly tuning the right parameter. We notice that KPCA provides good flexibility thanks to the choice of any possible kernel which fits the data geometry. The most simple algorithms are the PCA, and the distance function MorphPCA. Then, we have the scale-space decomposition Mor-phPCA, and finally the pattern spectrum MorphPCA. The second issue analysed is the computational complexity, and the third one is the memory requirements. From a computational viewpoint, the most demanding step in the PCA is the SVD, which can be done in O(D 3 ). PCA is the technique with the smallest computational need. On the contrary, the computational requirement of KPCA is O(n 3 ); since n D, this kind of algorithm seems infeasible in standard hyperspectral images. That is the reason why most of hyperspectral KPCA techniques use tricks to be able to deal with the high number of spectra [START_REF] Mathieu | Kernel principal component analysis for the classification of hyperspectral remote sensing data over urban areas[END_REF][START_REF] Mohan | Spatially coherent nonlinear dimensionality reduction and segmentation of hyperspectral images[END_REF][START_REF] He | Using diffusion geometric coordinates for hyperspectral imagery representation[END_REF]. All these techniques lead to a spatial distortion, which is not avoidable by the need of a sampling procedure aiming at reducing the number of spectra. Between the complexity of PCA and KPCA, we have our proposed MorphPCA algorithms. Regarding MorphPCA, the computationally demanding step is the computation of the morphological representation used in the corresponding covariance matrix. The complexity estimation has been carried out each time in the worse case, however efficient morphological algorithms can improve this part. Distance function MorphPCA is last demanding, then the scale-space decomposition MorphPCA and finally the pattern spectrum MorphPCA. Regarding memory needs, for the PCA, the pattern spectrum MorphPCA, and the CHAPTER 2. MORPHOLOGICAL PRINCIPAL COMPONENT ANALYSIS FOR HYPERSPECTRAL IMAGE ANALYSIS 63 distance function MorphPCA, the steps requiring more memory is the storage of the covariance matrix, just of O(D 2 ). The Spatial/Spectral MorphPCA needs to store 2 covariance matrices then its memory need is O(2D 2 ); similarly the scale-space MorphPCA needs to store 2S + 1 covariance matrices, then its required memory is O((2S + 1)D 2 ). Note that KPCA uses a Gram matrix of size (n × n).

Technique Parameter Computational Memory (1) (2) (3) PCA Prop O(D 3 ) O(D 2 ) MorphPCA Morpho-1 Prop, S O(DnS S (2S + 1)) O(D 2 (2S + 1)) MorphPCA Morpho-2 Prop, S O(DnS S (2S + 1)) O(D 2 ) MorphPCA Morpho-3 Prop O(Dn(b -a)) O(D 2 ) MorphPCA Morpho-4 β Prop, β O(DnS S (2S + 1)) O(2D 2 ) KPCA Prop, K O(n 3 ) O(n 2 )
Table 2.1: Comparison of the properties of dimensionality reduction algorithms for hyperspectral images.

Evaluation on hyperspectral images

The assessment of the performance of PCA and MorphPCA has been carried out on three hyperspectral images. The first image was acquired over the city of Pavia (Italy) and it represents the university campus. The dimensions of the image are 610 × 340 pixels, with D = 103 spectral bands and its geometrical resolution is of 1.3 m. We also used a second hyperspectral image which represents the University of Houston campus and the neighbouring urban area at the spatial resolution of 2.5 m and which dimensions are 349 × 1905 pixels and D = 144 spectral bands [START_REF] Debes | Hyperspectral and lidar data fusion: Outcome of the 2013 grss data fusion contest[END_REF].

The third image, acquired over the region of the Indian Pines test site in Northwestern Indiana, is composed for two-thirds of agriculture, and one-third of forest. The dimensions of this image are 145 × 145 pixels, D = 224 spectral bands and its geometrical resolution is of 3.7 m. We have applied classical PCA and the different variants of MorphPCA to Pavia hyperspectral image. Figure 2.12 shows the first three eigenimages, visualized as a RGB false color. We note that the pattern spectrum MorphPCA requires d = 5 to represent 92% of the variance whereas the other approaches only impose d = 3. An interesting aspect observed on the projection of the 103 spectral channels of Pavia hyperspectral image into the first two eigenvectors is how PCA and the scale-space decomposition MorphPCA cluster the bands linearly, see Figure 2.13(a) and 2.13(b). Bands close in the projection are also near in the spectral domain, whereas the pattern spectrum MorphPCA, 2.13(c), and distance function MorphPCA, 2.13(d), tend to cluster spectral bands which are not necessary spectrally contiguous. This can be explained thanks to Figure 2.10, where the MorphPCA correlation matrices are different from the classical PCA one.

It can be noticed that in classical manifold learning techniques, the goal is to decrease the dimension of the data while keeping some properties on the data manifold. We work here on the manifold of the channels. This manifold is easier to use, but finding the good β in V Morpho-4 β that would maintain some properties of 2.4. MORPHPCA APPLIED TO HYPERSPECTRAL IMAGES the manifold is not always easy, since we had to deal with a double optimization problem, i.e., β and d. From a quantitative viewpoint, one can see in Table 3.1 that globally Morph-PCA produces a more homogenous regularization of the image than classical PCA, especially the distance function MorphPCA and Spatial/Spectral MorphPCA with an appropriate β = 0.2, which gives the lowest values of Error Homg . We noted that Error sparse spatially follows a different ranking. A good method is the one with a good trade-off between both criteria, since one wants a DR to be trustworthy, which is evaluated thanks to Error Homg . But if the signal is too noisy, one may prefer a sparser representation. According to these criteria, the distance function Morph-PCA and the pattern spectrum MorphPCA seem to have the best result. We can also note that if we use manifold learning parameters for criterion C3, see Figure 2.14, the pattern spectrum MorphPCA has the best results. With respect to criterion C5, we have computed the 3D pattern spectrum distribution of Pavia hyperspectral image and of the different reduced images into d components, see Figure 2.15. From this result, we can see that both PCA and scalespace decomposition MorphPCA follow very well the hyperspectral image, since their spatial and spectral cumulative distributions are similar. However if one would like to denoise the hyperspectral image thanks to a DR technique, these results are not always positive. If we compare the spatial and spectral cumulative distribution of the distance function MorphPCA and the one of the hyperspectral image, we notice that for small 3D size (i.e., small spatial/spectral variations) that can be considered as noise, there are differences between these two distributions. But when the size increases, the distribution of the the distance function MorphPCA tends to the hyperspectral image one. So it seems that the distance function MorphPCA simplifies the spectral/spatial noise, considered as the small 3D objects, but keeps the objects of interest. Finally, Table 3.2 summarizes the results of supervised classification of respectively Pavia and Indian Pine hyperspectral images. We note that results for Pavia image are quite similar in all the cases even if MorphPCA seems to be better than PCA. Therefore, we have chosen to focus on the Indian Pine image, which is more challenging for supervised classification benchmark, see also Figure 2.16 and 2.17. We note that MorphPCA improves the results, especially the scale-space decomposition MorphPCA. To evaluate the classification results, first we fixed the dimension d of the reduce image. We have chosen d = 5. Then, we used the least square SVM, which is a multi-class classification technique, contrary to classical two-class SVM.

V V Morpho-1 V Morpho-2 V Morpho-
V V Morpho-1 V Morpho-2 V Morpho-
We also used two simple kernels: the linear one, which is the simplest one, and the rbf one, which is appropriate for hyperspectral images since we can assume that these data follow Gaussian distribution. Finally, we study the influence of dimension d and of the size of the training set on the different DR techniques over the classification results. For this purpose we have depicted in Figure 2.18 and Figure 2.19 the evolution of the kappa statistics. From the latter plot we can see that the PCA and the pattern spectrum PCA have the worst results. By combining the spectral and the spatial information, a better classification can be achieved. This is the case of the distance function MorphPCA, the scale-space decomposition MorphPCA, and the Spatial-Spectral MorphPCA for β = 0.2. 2.5. CONCLUSIONS

Conclusions

We have shown in this chapter how to introduce spatial information in the process of dimensionality reduction thanks to mathematical morphology operators. The representation techniques that we introduced in the chapter are based on the notion of the MorphPCA. As we said in the introduction, it might be possible to consider for such purpose a Kernel Principal Component Analysis (KPCA), where the kernel handles jointly spatial and spectral information. However as discussed in Section 4.2, KPCA needs a Gram matrix of size n × n, where n is the number of pixels. It is impossible to manipulate such a matrix with our test images on a standard computer.

That is why there are some works trying to approximate the gram matrix in order to perform the KPCA. Some of them consider a small subset of the data which can be chosen randomly or according to some spatial information. This kind of approach has been used with hyperspectral images for spectral-spatial processing. But these kind of works approximate the kernel, to perform dimensionality reduction on the considered Hilbert space. Here dimensionality reduction is done on the space of the data without any approximation on the morphological covariance V Morpho .

That is why we proposed techniques that are simple in computation and in memory storage and that can reduce the dimension while considering the spatial information. To assess these techniques we used typical criteria of dimensionality reduction which evaluate the fact that some properties are kept on the manifold of data after the dimensionality reduction. Moreover we also proposed some criteria to evaluate the quality of the image after the dimensionality reduction. Some of them are based on mathematical morphology, namely the 3D pattern spectrum and the α-flat zone to check that the reconstructed image preserves global and local similarity to the original hyperspectral image. Finally, we also perform a classification of the reduced data with different techniques. According to the entire set of criteria, adding spatial information improves the dimensionality reduction. However as we can see a good dimensionality reduction is obtained when we combine spatial and spectral feature space. A technique that seems to fulfil this optimum is the MorphPCA based on the distance function.

Finally, PCA has multiple applications on image processing; typically one can use the PCA to perform denoising. For example in the case of multiple images representing the same scene but corrupted by a Gaussian noise (like on multi-temporal images), it is possible to use the PCA to reduce the dimension of the temporal data and then to project the data back on the high dimensional space, so that we reduce the influence of the noise. This process can be done with MorphPCA. Moreover another case would be to use MorphPCA on fusion of information techniques like pansharpening, which consists of increasing the spatial resolution of a multispectral or hyperspectral image thanks to a grey scale image at high resolution. Some techniques for pansharpening are based on PCA. In summary, MorphPCA can be an appropiate alternative to PCA in different image processing applications.

INTRODUCTION

Introduction

The supervised classification of pixels in images is a difficult task, that consists in associating to each observation, which are pixels, a class that is not observed. This problem is common in remote sensing imaging, here we will apply it on multi/hyperspectral images. Hyperspectral images consist of very high-dimensional pixel observations that allows reconstruction of the spectral profiles of objects imaged thanks to the acquisition of several hundred narrow spectral bands. The supervised classification of these pixels is a challenging task, which commonly arises in remote sensing imaging [START_REF] Camps-Valls | Advances in hyperspectral image classification: Earth monitoring with statistical learning methods[END_REF][START_REF] Gualtieri | Support vector machines for hyperspectral remote sensing classification[END_REF][START_REF] Fauvel | Spectral and spatial classification of hyperspectral data using svms and morphological profiles[END_REF][START_REF] Fauvel | Advances in spectral-spatial classification of hyperspectral images[END_REF], but also in other domains. The structure of hyperspectral imagery imposed at first a lot of research to focus on spectral information without considering the spatial arrangement of pixels on a regular grid. The use of spatial information is a useful to fight the curse of dimensionality explained [13], and so the spectral variability of the vectors. This variability reduces the classification performance, and makes the distances between spectra meaningless. Moreover, neighboring samples might not be independent. So that if a given pixel belongs to a class, its neighbors might likely have the same class. These kinds of assumptions gave birth to the concept of regionalized random variable [START_REF] Matheron | Traité de géostatistique appliquée[END_REF] which is a key notion of geostatistics. To be able to take into account this spatial information, we focus here on the local texture of the image. Hence, we propose a novel approach to classification able to add this spatial information in order to describe the local texture. Our technique is based on kernel embeddings of distributions and on the scattering transform and they both use both the spatial and spectral information of the data. First we present the data and we introduce the scattering transform in Section 3.3, then we provide in Section 3.4 the background on kernel embeddings of distributions, random features for fast approximations to kernel methods. In Section 3.6 we describe the kernel mean map and study the consistency and convergence rate of the proposed method and empirical evaluation is given in Section 3.9. We note that this chapter is an extended and improved version of our following contributions [START_REF] Franchi | Hyperspectral image classification with support vector machines on kernel distribution embeddings[END_REF][START_REF] Franchi | A deep spatial/spectral descriptor of hyperspectral texture using scattering transform[END_REF].

Related works

Many techniques aim to include the spatial information in the image classification process. Of particular interest are those combining feature space representations describing the spatial information with those describing the pixels. Morphological feature spaces have been considered in several publications, with impressive results [START_REF] Dalla Mura | Morphological attribute profiles for the analysis of very high resolution images[END_REF][START_REF] Fauvel | Spectral and spatial classification of hyperspectral data using svms and morphological profiles[END_REF][START_REF] Pesaresi | A new approach for the morphological segmentation of high-resolution satellite imagery[END_REF][START_REF] Lefevre | Hyperspectral image classification from multiscale description with constrained connectivity and metric learning[END_REF]. On the other hand, kernel methods have also been studied extensively, and more particularly the compositions of kernels [START_REF] Li | Generalized composite kernel framework for hyperspectral image classification[END_REF][START_REF] Li | Multiple feature learning for hyperspectral image classification[END_REF][START_REF] Camps-Valls | Composite kernels for hyperspectral image classification[END_REF][START_REF] Mathieu | Kernel principal component analysis for the classification of hyperspectral remote sensing data over urban areas[END_REF], which allow building new feature space representations. We marry these approaches with a framework of [START_REF] Smola | A Hilbert space embedding for distributions[END_REF][START_REF] Muandet | Learning from distributions via support measure machines[END_REF][START_REF] Szabó | Two-stage sampled learning theory on distributions[END_REF][START_REF] Oliva | Deep mean maps[END_REF], where instead of the usual feature map, sending each data point to the feature space, a whole distribution can be represented in the RKHS. This idea yields a framework for learning on distributions via their representations in this RKHS. In our approach, each pixel is associated to an innovative image representation technique called the scattering transform [START_REF] Mallat | Group invariant scattering[END_REF][START_REF] Bruna | Invariant scattering convolution networks[END_REF]. It has successfully been used to classify images [START_REF] Sifre | Rotation, scaling and deformation invariant scattering for texture discrimination[END_REF] and sounds [START_REF] Andén | Multiscale scattering for audio classification[END_REF], and in this work we adapt it to be able to perform spatial pixel-wise classification. So on our approach, each pixel is described as a distribution of its neighbours of scattering transform CHAPTER 3. INVARIANT SPATIAL CLASSIFICATION OF MULTI/HYPER-SPECTRAL IMAGES 75 embedded on another RKHS. Another related line of work is that of [START_REF] Volpi | Spatially aware supervised nonlinear dimensionality reduction for hyperspectral data[END_REF], where the mean map is used on hyperspectral data to perform a dimensionality reduction. There are also recent works based on deep learning patterns [START_REF] Chen | Selected Topics in Applied Earth Observations and Remote Sensing[END_REF][START_REF] Ma | Hyperspectral image classification via contextual deep learning[END_REF][START_REF] Liang | Hyperspectral imagery classification using sparse representations of convolutional neural network features[END_REF], with interesting results, however they cannot compete yet against morphological features. Our proposal of deep learning technique based on texture description has the same order of classification accuracy than morphological descriptors. Dans ce chapitre, une nouvelle approche pour la classification des pixels d'images multi / hyper-spectral en s'appuyant à la fois sur l'information spatiale et spectrale des données. La méthode introduite repose sur un ensemble de technique récemment proposé pour pour apprendre la distribution des descripteurs de texture. Ces descripteurs visent à représenter le contenu de l'image tout en considérant les invariances liées à la texture et à ses transformations géométriques. En outre, nous considérons des invariances spectrales qui sont liées à la composition physique des pixels.

Invariant classification on hyperspectral images

Notation

Let E be a subset of the discrete space Z 2 , which represents the support space of a 2D image and F ⊆ R D be a set of pixels values in dimension D. Hence, it is assumed in our case that the value of a pixel x ∈ E is represented by a vector v ∈ F of dimension D. This vector v represents the spectrum at the position x, and D is the number of channels. Thus we write f : E → F the function associated to the image, such that f (x) = v.

Invariance properties of hyperspectral data

To improve the process of classification, the data should be invariant under one or more transformations underlying the "data" acquisition. For example when one want to classify the MNIST data set [START_REF] Lecun | The mnist database of handwritten digits[END_REF], the different letters of the data set are rotated, translated, and noised. So one would like to assign the class of images independently of the position or orientation of the characters. Such transformations may degrade the process of classification. That is why one would like to learn the invariant features of each individual class. To learn these invariants, different techniques can be used.

• One can increase the training set using variations of the training pattern according to all possible transformations. If the new training set that is synthesized is conform to the data, then the process of classification is improved.

• One can find a descriptor that may extract the invariances of each class. By extracting this feature space, the process of classification is more robust.

In this work we consider both approaches. We will first consider the second approach applied on the case of hyperspectral texture. The notion of texture [START_REF] Haralick | Textural features for image classification[END_REF] is linked with how human vision can interpret and recognize image structures, in particular, how human visual system isolate and segment textures. Often natural textures keep a pattern that periodically repeat itself on a certain scale. These pixels zones have a kind of statistical periodicity. Another very important property of textures is that they generally have multiple scales. Thus the size of the neighbourhood is an important information when dealing with texture information. For our part we consider that a texture is a set of connected pixels that appears as a homogeneous area. So these pixels follow the same distribution. In particular we work with hyperspectral textures. Then, the question is how can we learn an invariant descriptor of these hyperspectral textures. We consider that image f : E → F is a random function. To classify a pixel f (x) = v we will not only consider the pixel information, but also a patch around the position x, see Figure 3.1. Let us write P (s) v (f ) a patch surrounding the pixel v, of scale s. Then we would like to classify patches, to be able to attribute the class to the central pixel v. Because of this way of classifying, there would be huge redundancy, and we would like to learn invariants to describe these patches in this context.

In our case it is necessary for hyperspectral textures that characterize pixels to be invariant to changes in illumination conditions. Therefore we would like for our descriptor to be invariant to different kind of transformations.

Spatial invariance

Let us consider two patches, {P v1 (f ), P v2 (f )}, which can be seen as two images. Since we are first interested in the spatial invariance, we will assume that these patches are two grey scale images. These images are represented by mapping ψ(P v1 (f )) ∈ H, where ψ(P v1 (f )) is a texture descriptor of the patch, and H is a Hilbert space that characterize texture information. As explained on [START_REF] Bruna | Invariant scattering convolution networks[END_REF][START_REF] Mallat | Group invariant scattering[END_REF], we would like the texture descriptor to be invariant to rotation, translation, and to be Lipschitz continuous to deformation T , i.e. :

ψ(P v1 (f )) -ψ(P v2 (f )) 2 ≤ C f 2 sup x (T (x)) 2 , with P v2 (f ) = T (P v1 (f )), (3.1)
where C is a constant, and (T (t)) is the deformation gradient, such that its norm 2 expresses the magnitude of the deformation. This formula expresses a stability to CHAPTER 3. INVARIANT SPATIAL CLASSIFICATION OF MULTI/HYPER-SPECTRAL IMAGES 77 deformation.

Spectral invariance

To deal with spectral invariance we use a linear model of hyperspectral images. Pixel values of a hyperspectral image can be seen as embedded into a low dimensional set in a set of materials present in the scene, that are called endmembers. The dimension of each pixel can be reduced to the number of materials, where the new barycentric coordinates of each pixel are the abundances. By abundance we mean the positive quantity of each material on the pixel. Under the linear model, each pixel of the image can be written as a nonnegative combination of the different endmembers.

Thus, if we consider the spectrum at a pixel as the vector f (x) ∈ R D , then it can be written as:

f (x) = R r=1 a r (x)m r + n x , (3.2)
where {m r } R r=1 , m r ∈ R D represents the set of R endmembers, a r (x) the abundance at vector i of each endmember r, and n i an additive noise. This last term can be neglected. Hence, the extraction of endmembers can be seen as finding the simplicial cone containing the data. In general, for a given set of vectors there are many possible simplicial cones containing the vectors [START_REF] Donoho | When does non-negative matrix factorization give a correct decomposition into parts? In Advances in neural information processing systems[END_REF]. A way of reducing the number of possible representations consists in restricting the nonnegative coefficients a r to be a convex combination such that R r=1 a r (x) = 1, ∀x. By projecting all the data into the simplex of the abundance, each patch is then projected on the sample simplex. Then we have for each pixel v

P v (f ) ∈ R n 1 ×n 2 ×D -→ P v (a) ∈ R n 1 ×n 2 ×R (3.3)
With n 1 × n 2 being the spatial size of the patches. If we consider a second patch such that: P v2 (f ) = α.P v1 (f ), where α corresponds to an illumination change factor. The additional convex constraint guarantees that we are invariant to illumination changes, since P v2 (a) = P v1 (a)

Support vector machine classification for remote sensing

The classification of hyperspectral pixels turns out to be difficult, due to the high dimensionality of the pixels. In addition, the spectral data are noised, thus not all the sources that are present on the sample are necessary present on the acquired image. The process of image acquisition may also produce nonlinearities in the data. For all that reasons, it is considered that the data lie in a low dimensional complex manifold. The data present a high spatial and spectral variability. We represent in Figure 3. Let us consider a set of data {(v 1 , y 1 ), . . . , (v n , y n )} ∈ F × Y independent and identically distributed (i.i.d.) random variables following an unknown joint distribution P. We call this collection of points the training set. To simplify the discussion, we will consider a two-classes y i ∈ {-1, 1} situation, the case of multiple classes can be adapted thanks to well known techniques [START_REF] Bishop | Pattern recognition[END_REF]. Then, we consider that v i is a hyperspectral spectrum, and y i is equal to one of the classes. Let us write g: X -→ R a prediction function, such that if g(v i ) ≥ 0 then the class of v i is ŷi = 1, and in the contrary if g(v i ) < 0 the class of v i is ŷi = -1. The goal of the classification process is to find g such that g(v i )y i ≥ 0 for most of the data of the training set. Then the goal of the process is to maximize the margin function g(v i )y i , which is the smallest distance between the decision hyperplane and any data on the training set. Trying to find the hyperplane that would fit all the data is not always optimal, so instead of using this boundary, usually we determine a hyperplane by considering that only the nearby data matters. This particular points are called the support CHAPTER 3. INVARIANT SPATIAL CLASSIFICATION OF MULTI/HYPER-SPECTRAL IMAGES 79 vectors. Then the new margins to be maximized are those of the closest data v i from the possible hyperplane. This process of classification is called the Support Vector Machine (SVM). There are different kind of processes where the relevance of the training data set is not the same. This importance is modelled by a non increasing risk function Φ, which is a Lipschitz function. The choice of this risk function changes the algorithm, for example for the SVM the risk function is equal to Φ(u) = max(1u, 0). The objective function we want to minimize is :

R(g) = inf g∈H E (v,y)∼P (Φ(g(v).y)) . (3.4)
Since P is unknown the real function we optimize is the following empirical cost function :

R(g) = inf f ∈H 1 n . n i=1 (Φ(g(v i ).y i )) (3.5)
We assumed that the training set is linearly separable, so we can find an hyperplane to separate the two classes. Unfortunately this is not always the case, that is why there have been a lot of works in classification aiming at finding the good hyperplane that would bring the data into a space in which the classes are linearly separable. Moreover, in some cases, the classes may overlap so an exact separation of the classes would lead to a poor generalization. Then the SVM formulation is modified such that it allows for some points of the training set to be on the wrong side of the margin. This new constrain added to the SVM is controlled by a parameter C, see [START_REF] Friedman | The elements of statistical learning[END_REF] for more explanations. In the particular case of pixelwise classification of remote sensing images, an additional problem may appear: the size of the training set is small with respect to the variability of the data. To overcome this issue it is possible to generate data as proposed in [START_REF] Decoste | Training invariant support vector machines[END_REF].

Invariance thanks to training set generation

Learning an hyperplane to separate the data may be really tricky, and may depend on the number of available data. To encourage the generalization of the model, one can decide to generate random variation of the training set. The training set must be increased thanks to desired invariance. Then it is expected that thanks to this new training set the process will learn the invariances from the artificially enlarged training data, and will take into account them to build the hyperplan. The method of generating virtual examples must be applied carefully. As we discussed before, all the data do not have the same importance. Thus in [START_REF] Schölkopf | Incorporating invariances in support vector learning machines[END_REF], it has been decided to first proceed to a classical SVM, then learn the support vectors, and apply to this data the generation of random variation to increase the invariance. Finally another SVM is done with the new training set. This process is called the virtual SVM (VSVM). It has been applied on the MNIST data set [START_REF] Schölkopf | Incorporating invariances in support vector learning machines[END_REF], and also on remote sensing images [START_REF] Izquierdo-Verdiguier | Encoding invariances in remote sensing image classification with svm[END_REF]. On [START_REF] Izquierdo-Verdiguier | Encoding invariances in remote sensing image classification with svm[END_REF] the authors suggest some invariance properties for remote sensing. It happens that our feature space already does a patch-based classification, and is also invariant to scales. So we suggest to add Gaussian random noise to the SV, and also to increase the data set by applying some small rotation to the SV.

Invariance thanks to morphological profiles

Mathematical morphology operators are non-linear image processes based on the spatial structure of the image. Let f be a grey scale image which can be represented by a function. The two basic operators in morphology are the grey-level erosion and the grey-level dilatation whose definition are respectively given by [START_REF] Serra | Image Analysis and Mathematical Morphology[END_REF]:

ε b (f )(x) = inf h∈E {f (x -h) -b(h)}, (3.6) δ b (f )(x) = sup h∈E {f (x -h) + b(h)}, (3.7)
where b is a structuring function, which introduces the effect of the operators by the geometry of its support as well as the penalizations. We note that there are just convolutions in the (max, +)-algebra and its dual algebra [START_REF] Angulo | max, min)-convolution and mathematical morphology[END_REF]. We consider for simplicity uniform structuring functions which are formalised by their support set or shape B, called structuring element : b

(x) = 0 if x ∈ B -∞ otherwise
. By concatenation of these two basic morphological operators it is possible to obtain more evolved filters such as the opening and the closing [START_REF] Serra | Image Analysis and Mathematical Morphology[END_REF] :

γ B (f ) = δ B (ε B (f )) , (3.8) ϕ B (f ) = ε B (δ B (f )) .
(3.9)

These operators remove from f all the bright (opening) or dark (closing) structures where the structuring element B cannot fit. However they also modify the value of pixels where B fits. Thus to avoid these artefacts it has been proposed in [START_REF] Pesaresi | A new approach for the morphological segmentation of high-resolution satellite imagery[END_REF] to use geodesic opening and closing. Then by considering a set {γ

(i)
R }, i = 1...n, of indexed geodesic openings ,and a set {ϕ (i) R }, i = 1...n, of indexed geodesic closings where, typically, the index i is associated to the size of the structuring element. Then thanks to the granulometry axiomatic [START_REF] Serra | Image Analysis and Mathematical Morphology[END_REF] we obtain a scale-space representation of an image, which allows an image structure decomposition at different scales. Then the Morphological Profile (MP) of a grey scale image f at pixel x is defined as the 2 × n + 1 dimension vector:

M P (x) =             γ (1) R (f )(x) . . . γ (n) R (f )(x) f (x) ϕ (1) R (f )(x) . . . ϕ (n) R (f )(x)             (3.10)
To be able to use the MP on hyperspectral images, we first reduce the dimension of the data thanks to PCA, and then project the data on a d dimensional space which is of smaller dimension than the original space. So the hyperspectral image is represented by d grey scale eigen-images, then on each of these images we calculate the MP and we concatenate them. Hence, the spatial feature space is of dimension d × 2 × n + 1.

CHAPTER 3. INVARIANT SPATIAL CLASSIFICATION OF MULTI/HYPER-SPECTRAL IMAGES 81

Outline of the deep weighted mean map scattering representation

We will introduce here the outline of our process of classification. We will develop in the next sections each independent step. To do the classification of the pixels of the hyperspectral images the data is projected onto a feature space. To build this feature space, we first reduce the dimension of the hyperspectral image by unmixing, which involves the computation of endmembers, and then project them into a simplex of smaller dimension whose image representation is called abundance maps. Once the data lie on a small dimensional simplex, we apply the scattering transform to the abundance maps. Instead of extracting patches and applying the scattering transform to each patch, we apply the scattering transform and keep the size of the data. This process increases the speed of calculation of the feature field. Then we use this feature space to estimate a random rbf feature space of finite dimension see section 3.4.2 for details. Thanks to this explicit feature space a weighted mean map can be computed [START_REF] Flaxman | Who supported obama in 2012?: Ecological inference through distribution regression[END_REF]. This descriptor represents the layer 1, and which is the input for another step of random feature space computation and weighted mean pooling as input feature space the descriptor of layer 2. We can iterate it as many times as necessary, in our case of application due to the size of the training set there were no need to iterate it. Finally we concatenate the descriptors of the different layers, and apply a VSVM on it.

Scattering transform

Scattering transform has been recently studied [START_REF] Bruna | Invariant scattering convolution networks[END_REF][START_REF] Mallat | Group invariant scattering[END_REF] to describe image texture.

Here we deal with the notion of local texture on digital images, so we need a space and frequency representation for discrete signals. Moreover, since we apparently do not know what is the scale of the data, a multiscale analysis of the signal is naturally required. A tool to fulfil this requirement already exists and is based on the wavelets. Let ψ ∈ L 2 (R) be an orthonormal wavelet, so that : We can define the 1D wavelet transform as:

∞ -∞ | ψ(x)| 2 |x| 2 dx < ∞. ( 3
[W ψ f ] (a, b) = 1 |a| ∞ -∞ ψ x -b a f (x)dx, (3.12)
where the bar denotes the complex conjugate, a is the scaling, and b the time (in the case of a one dimensional signal). Wavelets are translation covariant, and so they follow the translation of the signal. If we consider two textures ν i and νi such that νi (x) = ν i (x + h) then we have :

[W ψ νi ] (a, b) = [W ψ ν i ] (a, b + h).
The authors of [START_REF] Bruna | Invariant scattering convolution networks[END_REF][START_REF] Mallat | Group invariant scattering[END_REF] explained that to build a translation invariant representation it is necessary to introduce a linear or non-linear operator M that commutes with translation.

The operator M introduced is the mean pooling, which corresponds to a mean filter in a fixed window. Nevertheless, applying this operator to wavelets might give a translation invariant result equal to zero because of (3.11). To overcome this issue, mean pooling is applied on the modulus of the wavelets. On images, we deal with two dimensional signals, so the wavelet should be invariant to translation in all the different directions. To achieve this goal, it is proposed in [START_REF] Bruna | Invariant scattering convolution networks[END_REF] to use a 2D discrete wavelet transform which is a representation of 2D data according to 4 variables: dilation, rotation, and position. Let

x ∈ E, if f (x) ∈ L 2 (R 2
) is square-integrable on E, then the 2D discrete wavelet transform is defined as :

[W ψ f ] (j, b, θ) = 2 -j ∞ -∞ ψ (r -θ .2 -j .(x -b))f (x)dx, (3.13)
where r θ is the 2D rotation matrix. Let us denote ψ(ω) the Fourier transform of the mother wavelet, and we write : ψ2 -j r -θ (ω) = ψ(2 -j .r -θ .ω) = ψλ (ω), where λ = 2 -j .r -θ is an index related to a pair of dilation/rotation parameters. Then, the wavelet transform of f is {f ψ λ (x)} λ , which is redundant and does not have orthogonality property. Since the modulus of the wavelet transform is used, we deal with U f (x) = {|f ψ λ (x)|} λ and then a mean pooling is applied. So finally one has Sf (x) = { f ψ λ (x) 1 } λ which is just the norm 1 of the wavelet coefficient, and it provides information on the sparsity of the wavelet. In our case, the mean pooling is a simple average filter with a square kernel of size s. Moreover, let us write

U [λ]f (x) = |f ψ λ (x)|.
Then, one issue is that all the multiscale variability of the signal is not handled by U f (x). A solution, proposed by [START_REF] Bruna | Invariant scattering convolution networks[END_REF][START_REF] Mallat | Group invariant scattering[END_REF] is, before calculating the mean pooling, to apply a wavelet transform on the set U f (x), and CHAPTER 3. INVARIANT SPATIAL CLASSIFICATION OF MULTI/HYPER-SPECTRAL IMAGES 83

to take their modulus and their mean pooling. They iterate this operation in order to recover all the information to be represented. Thanks to this technique, the variability of the texture is scattered into different paths p = (λ 1 , λ 2 , . . . , λ m ), where

U [p]f (x) = | . . . |f ψ λ 1 (x)| ψ λ 2 (x)| . . . ψ λm (x)|.
It was shown empirically in [START_REF] Bruna | Invariant scattering convolution networks[END_REF], that usually one can limit itself on the second order path. Then one define two representation vectors: U which represents the scattering transform without mean pooling, and S including the mean pooling:

U (x) =      f (x) |f ψ λ 1 |(x) ||f ψ λ 1 | ψ λ 2 |(x) . . .      , S(x) =      E(f (x)) E(|f ψ λ 1 |)(x) E(||f ψ λ 1 | ψ λ 2 |)(x) . . .     
It was shown in [START_REF] Mallat | Group invariant scattering[END_REF] that for appropriate wavelets, the scattering transform has the following properties:

• Contractive : S(f 1 ) -S(f 2 ) 2 ≤ f 1 -f 2 2 ;
• Preserve the norms : S(f 1 ) 2 = f 1 2 ;

• Stable to deformation :

S(f 1 ) -S(f 2 ) 2 ≤ C f 1 2 sup t (T (t)) 2 with f 1 = T (f 2 ).
In particular, for the Haar wavelet and with the appropriate scaling and orientation properties, the scattering coefficients are equivalent to the SIFT descriptor [START_REF] Lowe | Object recognition from local scale-invariant features[END_REF]. We can also notice that with the Haar wavelet, the scattering transform has some links with the variogram of order 1, which is a geostatistic descriptor of spatial information. We will focus on this interpretation in the Section 5.

Deep mean map

Mean map kernel

Let k : X × X → R be a positive definite kernel. By Moore-Aronszajn theorem [START_REF] Berlinet | Reproducing Kernel Hilbert Spaces in Probability and Statistics[END_REF], there is a unique RKHS H of real-valued functions on X where

g, k(•, v) H = g(v), for all g ∈ H, v ∈ X ,
implying that k(., .) corresponds to an inner product between features and, in particular, k(v, v ) = k(•, v), k(•, v ) H . This means that k(•, v) can be viewed as a feature of v ∈ X . For many typical choices of kernels k, the RKHS H is infinite-dimensional. Now, let X denote a random variable following a distribution P. The mean map or the kernel embedding [START_REF] Smola | A Hilbert space embedding for distributions[END_REF][START_REF] Sriperumbudur | Hilbert space embeddings and metrics on probability measures[END_REF] of P is defined as:

µ P := E X [k(•, X)] = X k(•, v) dP(v), (3.14)
where the expectation is over H. For characteristic kernels [START_REF] Sriperumbudur | Universality, characteristic kernels and RKHS embedding of measures[END_REF], which include Gaussian rbf, Matern family and many others, this embedding is injective on the space of all probability distributions (i.e., captures information on all moments, similarly to a characteristic function). Further, if we are given two random variables,

DEEP MEAN MAP

X following the distribution P, and Y following the distribution Q, the inner product between the corresponding embeddings is given as

µ P , µ Q H = E X,Y [k(X, Y )], (3.15)
which is sometimes referred to as a mean map kernel. For a random sample {v 1 , . . . , v n }, drawn i.i.d. from P, we can define the empirical mean map:

µ P = 1 n n i=1 k(•, v i ), (3.16)
and for random samples {v 1 , . . . , v n } from P and {v 1 , . . . , v m } from Q, we obtain the empirical mean map kernel:

µ P , µ Q H = 1 nm n i=1 m j=1 k(v i , v j ).
(3.17)

Random features for kernels

The computational and storage requirements for kernel methods on large datasets can be prohibitive in practice due to the need to compute and store the kernel matrix. If we consider a dataset of n D-dimensional observations, the storage requirements are O(n 2 ) and the calculation takes O(Dn 2 ) operations. A remedy developed by [START_REF] Rahimi | Random features for large-scale kernel machines[END_REF] is to approximate translation-invariant kernels in an unbiased way using a random feature representation. Namely, any translation-invariant positive definite kernel k, such that

∀(v, v ) ∈ X 2 , k(v, v ) = κ(v -v ) can be written as k(v, v ) = E ω∼Λ cos(ω v) cos(ω v ) + sin(ω v) sin(ω v ) ,
where ω ∈ R D follows some distribution Λ (spectral measure of the kernel). Thus, by sampling i.i.d. vectors ω 1 , . . . , ω N from Λ, we can approximate kernel k by k defined by:

k(v, v ) = 1 N N j=1 cos(ω j v) cos(ω j v ) sin(ω j v) sin(ω j v ) ,
so that the original feature map k(•, v), potentially living in an infinite-dimensional space, is approximated by an explicit 2N -dimensional feature vector:

Ẑ(x) = 1 N cos(ω 1 v), . . . , cos(ω N v), sin(ω 1 v), . . . , sin(ω N v) T . (3.18)
Thus, the mean map and the mean map kernel can be estimated using these finitedimensional representations. In this part of thesis, we will focus on Gaussian rbf kernels for which the spectral measure Λ is also Gaussian.

Random features mean map on hyperspectral images

Let us now turn our attention to a hyperspectral image f . Around each pixel location x i , we consider a square patch P (s)

x i (f ) of size s where we will treat the pixels as a random sample from a distribution P i specific to the location x i . Instead CHAPTER 3. INVARIANT SPATIAL CLASSIFICATION OF MULTI/HYPER-SPECTRAL IMAGES 85 of calculating the kernel between individual data points, we will calculate the kernel between these distributions. An empirical mean map kernel is thus given simply by:

K mm (x i , x j ) = μP i , μP j H (3.19) = 1 s 2 l 1 ∈Px i l 2 ∈Px j k(f (x l 1 ), f (x l 2 )) ≈ 1 s 2 l 1 ∈Px i l 2 ∈Px j Ẑ(f (x l 1 )) Ẑ(f (x l 2 )),
where f (x) denotes the measurement vector at location x and in the last line we employ a random feature approximation of k.

It should be noted that there maybe outliers in a patch, which can damage the estimation of the mean. Similarly to the work of [START_REF] Flaxman | Who supported obama in 2012?: Ecological inference through distribution regression[END_REF], we proposed to use a weighted mean map, where the weights depend on spatial/spectral information. The kernels we obtain, called convolutional kernels, have also been used in [START_REF] Mairal | Convolutional kernel networks[END_REF]. In contrast to [START_REF] Mairal | Convolutional kernel networks[END_REF], however, we will use random features expansions to explicitly represent the feature space.

More precisely, the convolutional kernel is defined as:

KCN (x i , x j ) = l1∈Px i l2∈Px j w ij (β, x l1 , x l2 )e -1 2σ 2 f (x l1 )-f (x l2 ) 2 2 , (3.20)
The value of the weight we selected is

w ij (β, x l1 , x l2 ) = e -1 2β 2 x l1 -x i 2 2 .e -1 2β 2 x j -x l2 2
2 . This formula can be represented by Figure 3

Weighted mean map for a spatial regularization

Let us consider a sample {f (x 1 ), . . . , f (x n )} ∈ P (s)

x I (f ) drawn i.i.d. from P I . The empirical mean map is then:

µ P I = 1 n n i=1 φ(f (x i )) (3.21)
While µ P I is an often used estimator of µ P I ∈ H, due to the fact that the Hilbert space H is high-and potentially infinite dimensional whereas the number of data 3.4. DEEP MEAN MAP on the patch is small, then Stein's phenomenon implies that it is inadmissible. To improve the performance of the estimator, the authors of [START_REF] Muandet | Kernel mean estimation via spectral filtering[END_REF] propose to work with a family of "shrinkage" estimators. Given that the empirical mean map is the solution of the following minimization problem:

µ P I = argmin g∈H 1 n n i=1 φ(f (x i )) -g 2 H , (3.22)
the shrinkage mean map estimator μλ solve the regularized minimization problem instead:

μλ = argmin g∈H 1 n n i=1 φ(f (x i )) -g 2 H + λ. g 2 H . (3.23)
Here, λ > 0 is a shrinkage parameter that balances the bias and the variance of the model. The choice of λ is explained in [START_REF] Muandet | Kernel mean shrinkage estimators[END_REF]. Finally the class of estimators proposed in [START_REF] Muandet | Kernel mean shrinkage estimators[END_REF] is :

μλ = n i=1 β(λ) i φ(f (x i )), with β(λ) = g λ (K)K1 n , (3.24)
with g λ (γ) = 1/(λ + γ) where γ is an eigenvalue of K. So we have:

μλ = (K + nλI) -1 K1 n .(φ(f (x 1 )), ..., φ(f (x n ))) t . (3.25)
This kind of shrinkage is interesting but it does not consider the spatial consistency that is really important in spatial data such as image. So to handle this problem, we propose to work with a spatial shrinkage mean map estimator μS P I that solve the regularized minimization problem:

μS P I = argmin g∈H φ(f (x)) -g 2 H dP I (x)dG β,I (x) , (3.26)
where dG β,I is a positive measure expressing the spatial relation between the data which are hyperspectral pixels. We decided to fix

dG β,I (x) = e -1 2β 2 x-x I 2 2 dx.
The parameter β is difficult to fix by a prior knowledge, that is why we learn it during the classification. Since we do not have the distribution P I , we work with the empirical mean map. Therefore, we have:

μS P I = argmin g∈H 1 n n i=1 φ(f (x i )) -g 2 H e -1 2β 2 x i -x I 2 2 . (3.27)
By deriving according to g, we obtain that

µ P I = 1 n i=1 e -1 2β 2 x i -x I 2 2 n i=1 φ(f (x i ))e -1 2β 2 x i -x I 2 2 . (3.28)
This kind of mean map is structurally similar to the bilateral filter, that approximate the non-linear adaptive diffusion process. CHAPTER 3. INVARIANT SPATIAL CLASSIFICATION OF MULTI/HYPER-SPECTRAL IMAGES 87

Geostatistics of the feature field

Let us consider a probability space (Ω, A, P ) and a domain D ∈ E. A random field or random function on the spatial domain D with values in F ∈ R d is a function of two variables, denoted Z(x, ω). In our case, we will simplify the random field by considering that d = 1, which is the case that we study. For each x 0 ∈ D, Z(x, :) : ω -→ Z(x, ω) is a random variable on (Ω, A, P ). Moreover, for each ω 0 ∈ Ω, Z(: ω 0 ) : x -→ Z(x, ω 0 ) is a function of D -→ F . We will write Z(x) the random function at the position x. To study the properties of this function, one may seek to characterize stationary properties. To measure the variability of the random function Z(x) at different scales, we can calculate a measure of dissimilarity between two positions x 1 and x 2 . This dissimilarity between two values, designated by γ * , is defined by:

γ * = 1 2 (Z(x 1 ) -Z(x 2 )) 2 (3.29)
To better characterize the fields, the dissimilarity γ * is made dependent on the distance between the two points and orientation. By averaging dissimilarities of γ * for all values between N H pairs of points connected by a vector h, we get the notion of experimental variogram of order 2:

γ * (x i , h) = 1 2N H N H i=1 Z(x i ) -Z(x i + h) 2 2 . (3.30)
Usually, it is observed that the similarity values increase on average depending on the spatial distance of the measurement points and frequently leveled variation at huge distances. Behavior at very small scales, near the origin of the variogram, is of critical importance because it is an indicator of the degree of continuity of the regionalized variable. Furthermore, in the case of an image where it is hardly conceivable to have stationarity and isotropy properties, it seems important to limit our study of the variogram to work with low values of h. Finally, we often look at the variogram of order 2 to study the properties of a random field. But it is also quite interesting to look at the order variagram alpha, that bring other information, that is defined by:

γ α (h) = 1 2N H N H i=1 |Z(x i ) -Z(x i + h)| α .
(3.31)

Let us come back to our framework, with scattering transform, and especially with a Haar wavelet. We can see that the first order of the scattering transform U , at the first scale, corresponds to γ * (h = 1). Then to be able to deal with the variagram of order 1 we just need to calculate the local mean S. Then when we consider scattering coefficients at higher scales in a certain way we always focus on variagram of order 1 with h = 1, but instead of doing it on the original image, we do it on an image at a lower resolution. Then, since we decrease the resolution, each point on the random field at lower resolution is a combination of other points. So we "calculate" variogram of order 1 with more than two points. This can be seen as increasing the robustness of the calculation of the variogram, and is equivalent to increasing the 3.5. GEOSTATISTICS OF THE FEATURE FIELD length of h. Then let us consider that U (x) represents the concatenation of γ * (x, h) with different orientations and lengths of h.

Let us consider that we want to characterize the distribution near two positions x and y. Let us write P (s)

x (U ), P (s) y (U ) the two patches extracted at two positions with P (s)

x (U ) ∼ P x , and P (s) y (U ) ∼ P y . We can try to estimate the mean, the variance or other moments of P x and P y and compare them. Or we could try to use other estimator such as the widely used Kullback-Leibler divergence which would require the density estimation. However, since apparently we have no prior information on these distributions, we do not know if it would be sufficient. Many distances on distributions can be used. Here we focus on a innovative tool called the maximum mean discrepancy (MMD) which is defined as the squared distance between their embeddings in the RKHS [START_REF] Smola | A Hilbert space embedding for distributions[END_REF][START_REF] Song | Learning via hilbert space embedding of distributions[END_REF]

MMD(P x , P y ) = µ Px -µ Py 2 H . (3.32)
This distance is equivalent to finding the function on the RKHS that maximizes the difference in expectations between the two probability distributions, i.e.:

MMD(P x , P y ) = sup ||f || H ≤1 E X [f (P (s) x (U ))] -E Y [f (P (s) y (U ))] (3.33)
However, most of the time we work with a sampling of the distribution. Let us consider a random set of samples {v 1 , . . . , v n } from P x and {v 1 , . . . , v m } from P y . Then we can approximate the MMD by the empirical estimate of the MMD, defined by :

MMD(P x , P y ) = 1 n n i=1 φ(v i ) - 1 m m i=1 φ(v i ) 2 H (3.34) MMD(P x , P y ) = 1 nm n i=1 m j=1 k(v i , v j ) + k(v i , v j ) -2k(v i , v j ) (3.35)
where n and m are respectively the size of the patches P x and P y . Let us now study the link of the MMD with characteristic functions (Φ Px , Φ Py ) respectively of (P x , P y ), such that:

MMD(P x , P y ) = µ Px -µ Py 2 H MMD(P x , P y ) = k(v x , v x ) + k(v y , v y ) -2k(v x , v y ) dP x (v x )dP x (v x )dP y (v y )dP y (v y )
Since we work with positive definite and continuous kernels, we can use the Bochner theorem and rewrite the MMD as:

MMD(P x , P y ) = . . . e jω t (vx-v x ) + e jω t (vy-v y ) -2e jω t (vx-vy) Λ(ω)dP x (v x )dP x (v x )dP y (v y )dP y (v y )dω As we can see in Figure 3.7, thanks to the MMD, we calculate the difference between two characteristic functions of the variogram that have been filtered by a low pass filter Λ(ω). In our case the low pass filter is a centred Gaussian function of parameter σ, that is learned during the validation step. Thanks to this technique we expect that we are able to learn well the representation of the texture according to the class we want to discriminate.

Support vector machines on kernel distribution embeddings

In images, it is natural to assume that pixels close between them are more linked than pixels fare away. Then studying region properties to find the class of pixels 3.6. SUPPORT VECTOR MACHINES ON KERNEL DISTRIBUTION EMBEDDINGS seems a promising approach. We can furthermore consider that the image is a set of classes, and also that pixels are gathered in set of bags of pixels that follow the same distribution. These bags are composed of pixels close to each other. In [START_REF] Szabó | Two-stage sampled learning theory on distributions[END_REF] the author proposed to work with a super-pixels strategy to attribute to each super pixel the correct class. The difficulty with this kind of approach is that if the super-pixel assigned has a wrong class, then all the pixels on this super-pixel are wrongly classified. Therefore this technique depends not only on the feature extracted but also highly on how we select the super-pixels. Here we propose to use a sliding window approach, which has the advantage of a less important prior. We have calculated in [START_REF] Franchi | Hyperspectral image classification with support vector machines on kernel distribution embeddings[END_REF] that if we consider a sliding window approach or a super-pixel approach and we calculate a kernel mean map SVM, then the SVM converges to a SVM between the distributions of the super-pixels.

Let us consider that the data are partitioned into sets following the same distribution, then the structure of our data is given by {({v i,n

} N i n=1 , y i )} l i=1 with v i,1 , . . . , v i,N i i.i.d.
∼ v i , where (v i , y i ) are drawn from a joint meta distribution M. We follow the notation of [START_REF] Szabó | Two-stage sampled learning theory on distributions[END_REF]. We represent this problem in Figure 3

R(f ) = inf f ∈H E (v,y)∼M (Φ(f (v)y)) . (3.37)
where Φ is a loss function. We can modify it to mean map embedding classification problem:

R µ (f ) = inf f ∈H E (v,y)∼M (Φ(f (µ v )y)) . (3.38)
We can also write the empirical risk function, for mean map embedding classification problem:

Rµ (f ) = inf f ∈H 1 n n i=1 (Φ(f (µ v i )y i )) .
(3.39)

CHAPTER 3. INVARIANT SPATIAL CLASSIFICATION OF MULTI/HYPER-SPECTRAL IMAGES 91

Finally we can also write the empirical risk function, for the empirical mean map embedding classification problem:

Rμ (f ) = inf f ∈H 1 n n i=1 (Φ(f (μ v i )y i )) (3.40)
Then, we would like to obtain an inequality between R µ (f ) and Rμ (f ). To do that, inspired by [START_REF] Muandet | Learning from distributions via support measure machines[END_REF], we derive a inequality between R µ (f ) and R(f ).

Theorem 17 Given that x ∼ P an arbitrary probability distribution with variance σ 2 , a Lipschitz continuous function f : R → R with constant C f , an arbitrary loss function Φ : R → R that is Lipschitz continuous with constant C l , it follows that :

R µ (f ) -R(f ) ≤ C l C 2 f E (v) v -µ v 2 E (y) (y 2 ) (3.41)
Proof. First we have:

R µ (f ) -R(f ) ≤ E (v,y)∼M [Φ(f (v).y) -Φ(f (µ v ).y)] R µ (f ) -R(f ) ≤ E (v,y)∼M |Φ(f (v).y) -Φ(f (µ v ).y)|. Since Φ is Lipschitz continuous we obtain: R µ (f ) -R(f ) ≤ C l E (v,y)∼M |f (v).y -f (µ v ).y|.
Thanks to the Cauchy-Schwarz inequality, we finally get:

R µ (f ) -R(f ) ≤ C l E (v) (f (v) -f (µ v )) 2 E (y) (y 2 ) Since f is Lipschitz continuous we have: R µ (f ) -R(f ) ≤ C l C 2 f E (v) v -µ v 2 E (y) (y 2 )
Then we might use [30] where we have an inequality between R µ (f ) and Rµ (f ) :

Theorem 18 Let G = Φ(H, .) denote the loss class, let R n (G) denote the Rademacher complexity. Let Σ(G) 2 = sup g∈G E(g 2 ) be a bound on the variance of the functions in G. If the trace of the kernel is bounded, the loss function Φ : R → R that is Lipschitz continuous, for any δ > 0 , the following bound holds with probability at least 1δ

Rµ (f ) -R µ (f ) ≤ 8R n (G) + Σ(G) 8 log(2/δ) n + 3 log(2/δ) n
Theorem 19 Given that f : R → R is a Lipschitz continuous function with constant C f , an arbitrary loss function Φ : R → R that is Lipschitz continuous with constant C l , it follows that :

Rμ (f ) -Rµ (f ) ≤ 1 n C l C 2 f Ê µ v -μv 2 Ê (y) 2 Proof. Rμ (f ) -Rµ (f ) ≤ 1/n. [ n i=1 Φ(f (μ v i ).y i ) -Φ(f (µ v i ).y i )] Rμ (f ) -Rµ (f ) ≤ 1/n. ( n i=1 |Φ(f (μ v i ).y i ) -Φ(f (µ v i ).y i )|) Since Φ is Lipschitz con- tinuous we have: Rμ (f ) -Rµ (f ) ≤ 1/n.C l ( n i=1 |(f (μ v i ) -f (µ v i )).y i |).
Thanks to the Cauchy-Schwarz inequality we have:

Rμ (f ) -Rµ (f ) ≤ 1/n.C l ( n i=1 |y i | 2 ) ( n i=1 |f (μ v i ) -f (µ v i )| 2 ) Since f is Lipschitz continuous we obtain: Rμ (f ) -Rµ (f ) ≤ 1/n.C l C 2 f ( n i=1 (y i ) 2 ) ( n i=1 μv i -µ v i 2 )
We also need the following theorem proved in [START_REF] Smola | A Hilbert space embedding for distributions[END_REF]

KERNEL MEAN MAP SCATTERING AND INVARIANCE

Theorem 20 Assume that g ∞ ≤ R for all g ∈ H with g H ≤ 1, and that k is an universal kernel. Then with probability at least 1δ :

|µ[P ] -µ[X]| ≤ 2R n (H, P ) + R log (1/δ) /n
where R n (H, P ) denotes the Rademacher average associated with P and H.

Then by combining the previous results we easily have the following theorem.

Theorem 21 Given the conditions of the previous theorems, then with probability at least 1δ, we have:

Rμ (f ) -R(f ) ≤ C l C 2 f E (v) v -µ v 2 E (y) (y 2 ) + 2R n (H, P ) + R log (1/δ) /n Ê (y) 2 +8R n (G) + Σ(G) 8 log(2/δ) n + 3 log(2/δ) n
This theorem states that if the random variable v is concentrated around its mean and the functions f and Φ check the conditions of the previous theorems, then the loss deviation Rμ (f ) -R(f ) will be small.

Kernel mean map scattering and invariance

As we have discuss above, it has been established in [START_REF] Mallat | Group invariant scattering[END_REF] that for appropriate wavelets, the scattering transform has the following properties:

• it is contractive: S(f 1 ) -S(f 2 ) 2 ≤ f 1 -f 2 2 ,
• it preserves the norms: S(f 1 ) 2 = f 1 2 ,

• it is stable to deformations:

S(f 1 ) -S(f 2 ) 2 ≤ C f 1 2 sup t (T (t)) 2 with f 1 = T (f 2 ).
In our case, we perform the average pooling on the scattering coefficients U [p]f that have been projected on the rbf feature space. We denote this scattering transform S H [p](f ). Then one can wonder if these properties of invariance are still satisfied. We can prove that it is still contractive, and stable to deformation, but it does not preserve the norms anymore.

Proof. Since the exponential function is convex we have that:

e u ≥ u + 1 ∀u ∈ R Thus, e -v-v 2 2σ 2 ≥ - v -v 2 2σ 2 + 1 2 -2e -v-v 2 2σ 2 ≤ v -v 2 σ 2 CHAPTER 3. INVARIANT SPATIAL CLASSIFICATION OF MULTI/HYPER-SPECTRAL IMAGES 93
Then, replacing by the corresponding rbf feature space we have:

φ(v) -φ(v ) 2 =< φ(v) -φ(v ), φ(v) -φ(v ) > φ(v) -φ(v ) 2 =< φ(v), φ(v) > + < φ(v ), φ(v ) > -2 < φ(v), φ(v ) >= 2 -2e -v-v 2 2σ 2
.

Hence,

φ(v) -φ(v ) 2 ≤ v -v 2 σ 2
this means that the mapping of the rbf kernel is Lipschitz-continuous of parameter 1/σ 2 . Using this result we have:

S H (f 1 ) -S H (f 2 ) 2 = E (φ(U f 1 ) -φ(U f 2 )) 2
According to the fact that φ is Lipschitz-continuous, we obtain:

S H (f 1 ) -S H (f 2 ) 2 ≤ 1/σ 2 S(f 1 ) -S(f 2 ) 2
Finally one gets:

S H (f 1 ) -S H (f 2 ) 2 ≤ 1/σ 2 f 1 -f 2 2 (3.42)
and

S H (f 1 ) -S H (f 2 ) 2 ≤ C/σ 2 f 1 2 sup t (T (t)) 2 (3.43)
So the new descriptor is stable to deformation and contractive depending on a the parameter σ of the rbf kernel.

Multiple kernel mean map

Learning the hyperparameter of the rbf kernel is hard, and unfortunately the classification algorithm is really sensitive to this choice. That is why we first consider that [START_REF] Takeuchi | Nonparametric quantile estimation[END_REF]. However, it happens that this choice may not be perfectly adapted for

σ = median({ v i -v j R D , (i, j) ∈ [1, n] 2 }) as proposed in

EXPERIMENTS

Most of these techniques [START_REF] Rakotomamonjy | [END_REF] follow the large margin framework of the SVM: first they start with a set of training samples D = {(v 1 , y 1 ), . . . , (v n , y n )}. Then, they formulate their optimization problem by as:

min k∈K min f ∈H K λ f H K + n i=1 l(y i , f (v i )) , (3.44)
where l is a loss function, K is the optimization domain of the candidate kernels.

There are different ways to solve this optimization problem. This kind of technique lean on the fact that we can learn the coefficient β m thanks to a good training set. However here we will consider that we do not have enough information to learn the parameter thanks to the training set, and so we do not want to learn coefficients of the multiple kernel. We note that their are some links between the cost function of equation (3.44) and the cost function of the SVM on the sum kernels developed in [START_REF] Rakotomamonjy | [END_REF]. In any case, we are going to use the sum kernel. We use on the one hand the fact that the linear combination of kernels is still a kernel, such that on the other hand, we use universal kernels and the linear combination of universal kernels is still a universal kernel. Then we can perform multiple kernel mean maps. We note that the use of multiple kernel mean maps has already been proposed in [START_REF] Gretton | Optimal kernel choice for largescale two-sample tests[END_REF], but here we use it to classify data. Hence, we decide to first estimate σ thanks to the training and validation sets, we write this parameter σ . So we find the first kernel k 1 (v i , v j , σ ). Then we use as multiple kernel

K(v i , v j ) = k 1 (v i , v j , σ ) + k 1 (v i , v j , γ -1 σ ) + k 1 (v i , v j , γσ ),
where γ is a parameter that translate how much confident we are that the training set and the validation can generalize the testing set. Then we use this new kernel K to evaluate the mean map.

Experiments

Hyperspectral remote sensing

In this section, we evaluate the classification accuracy of the proposed approach using two hyperspectral images classically considered in the state-of-the-art: the AVIRIS Indian Pines data set, and the ROSIS University of Pavia data set. The first data set is an image of dimensions 145 × 145 pixels, with D = 224 spectral bands and its geometrical resolution is of 3.7 m. The dimensions of the second data set are 610 × 340 pixels, with D = 103 spectral bands and its geometrical resolution is of 1.3 m. On Pavia dataset there is a predefined testing set that we used in our experiments. In the first data set there is no common testing set so we generate 20 Monte-Carlo simulations, peaking randomly 5 pixels per class, then we aggregate the results of the supervised classification. We also generate training sets with 15 pixels per class, and with 50 pixels per class, however since there was not enough data for each class to have 50 points per class, we discarded the smallest classes. We use the C-SVM [START_REF] Chang | Libsvm: A library for support vector machines[END_REF] as classification algorithm. We learned the parameter C thanks to a validation on the training set. Then we compared our results to the the Morphological Profile (MP) [START_REF] Pesaresi | A new approach for the morphological segmentation of high-resolution satellite imagery[END_REF][START_REF] Dalla Mura | Morphological attribute profiles for the analysis of very high resolution images[END_REF] on pixel classification. We also use the product of the two kernels where one is the MP kernel and the other is the Kernel mean map (KMM kernel). This kind of technique has been previously explored in [START_REF] Li | Generalized composite kernel framework for hyperspectral image classification[END_REF][START_REF] Li | Multiple feature learning for hyperspectral image classification[END_REF][START_REF] Camps-Valls | Composite kernels for hyperspectral image classification[END_REF]. In contrast to previous works, which approximates the product of kernels thanks to addition of kernels, we can do the real multiplication since we work with finite dimension Hilbert spaces. Let us first illustrate the intermediary image images of our algorithm using Indian Pines. Figure 3.9 shows the scattering transform U before the weighted mean map. We also represent on Figure 3.9 the weighted mean map. As one may notice, both feature spaces provide spatial data that represents the structures present on the image.

From the results on Tables 3.1, 3.2 and, 3.3 we can note that the scattering transform can compete with morphological profiles [START_REF] Dalla Mura | Morphological attribute profiles for the analysis of very high resolution images[END_REF][START_REF] Li | Generalized composite kernel framework for hyperspectral image classification[END_REF] on Indian Pines data set. In addition, on the Pavia data set, we performed a Convolutional Neural Network (CNN) classification. CNN are mostly used to describe the texture of an image in order to classify it. Here we use CNN to classify patches. They are different CNN architectures that could be used. After trying several architectures we choose architecture inspired by the Oxfordnet [START_REF] Vedaldi | Matconvnet: Convolutional neural networks for matlab[END_REF]. The input patch is of size 32 × 32 × 3, such that to have a feature space of size 3, a PCA is performed on the Pavia hyperspectral image. We apply to these patches a first convolution of size 5 × 5 × 3 × 32. Then a max pooling with sub-sampling of factor 2, followed by a Relu layer is applied. Then we perform a new convolution layer of size 5 × 5 × 32 × 32 followed by a Relu and an average pooling with sub-sampling of factor 2. We apply a convolution layer of size 5 × 5 × 32 × 64, then a Relu layer and an average pooling with sub-sampling of factor 2 , followed by a new convolution layer of size 4 × 4 × 64 × 64. Then as regularizer we apply a dropout to cancel 30% to cancel the less powerful coefficients. Then we apply two fully connected layers. Obviously, this network is probably not the most optimal one, however it gives us a hint of some results on performance of convolutional neural networks applied on our problem of pixel classification. Moreover since the training set is quite small, the training set has been increased by rotating the patch according to 8 orientations : {0, π/4, π/2, 3π/4, π, 5π/4, 3π/2, 7π/4}. Table 3.4 gives the results of the different descriptors. One can conclude that our descriptor can compete with the others. From these two data set, we criticize the independence between the training set and the testing set since they are taken on the same image, these questions are studied in the following reference [START_REF] Liang | On the sampling strategy for evaluation of spectral-spatial methods in hyperspectral image classification[END_REF].

Multispectral remote sensing

We test now our machine learning on the Zurich dataset [START_REF] Volpi | Semantic segmentation of urban scenes by learning local class interactions[END_REF]. This data set was of 20 multispectral VHR images acquired over the city of Zurich (Switzerland), that has been studied in [START_REF] Volpi | Semantic segmentation of urban scenes by learning local class interactions[END_REF]. The typical image is composed of 1.2 million pixels and they are all composed of 4 channels RGB and near-infrared (NI). The spatial resolution is about 0.61 meters / pixel. The total number of classes is 8. There are no images having the 8 classes. On the contrary to [START_REF] Volpi | Semantic segmentation of urban scenes by learning local class interactions[END_REF], we work with pixels, instead of superpixels.

We did not apply any preprocessing on the multispectral images that is composed of R,G,B and NI channels. The multiclass SVM algorithm used corresponds to the liblinear library [START_REF] Fan | Liblinear: A library for large linear classification[END_REF] which is able to handle huge set of data with support vector classification proposed by Crammer and Singer [START_REF] Crammer | On the algorithmic implementation of multiclass kernel-based vector machines[END_REF]. However, because of the size of the dataset, we were not able to use the full potential of our descriptors. We just use the first order of the scattering transform with 3 rotations and 4 scales parameters, and also one layer of weighted mean map. Moreover the dimension of the estimated rbf feature space was just fixed to 30. This may degrade considerably the results of classifications.

To fix the different parameters we proceed to a validation. In order to have a training/ testing set independent, we first select 1 image that represents the testing

Conclusion

In this chapter we dealt with a hyper/multi-spectral pixel image classification in the context of remote sensing. We proposed a deep spatial and spectral descriptor that is able to learn the representation of the local texture of multivariate remote sensing images. Our deep descriptor has been assessed on 3 different data sets where the results were quite good. We evaluate the asymptotic properties of learning the local distribution thanks to this descriptor and a margin classifier. We compare this descriptor with different fields related to spatial statistical learning. This descriptor computes a translation invariant representation of the texture which is also Lipschitz stable to deformations. To evaluate this descriptor we just calculate the mean of the scattering coefficients provided by each layer embedded on the Hilbert space, of a rbf kernel. This links with mean map embedding. This classifier simplifies convolutional neural network since less parameters are needed. This can be an interesting way to 100 3.10. CONCLUSION perform deep learning when the training data set is limited. The proposed approach improves the classification of pixel, in the case of remote sensing. Moreover it could be interesting to see its results on other data sets. Another way to improve this technique might be to put more parameters when we perform the average weighted. This might be linked to the Perceptron neural network algorithm. We could also work on simulation of data. This might be really interesting and might be considered carefully.

Part III Fusion of information for multimodal SEM images

Chapter 4

Multimodal Scanning Electron Microscopy Images

Résumé Dans ce chapitre nous introduisons les concepts liés au projet effectué avec IFP Energies Nouvelles. Durant cet projet des images de microscope électronique à balayage (MEB) sous différentes modalités ont étudié. Ainsi ce chapitre vise à présenter ces différentes modalités.

Background on Scanning Electron Microscope (SEM)

The Scanning Electron Microscopy (SEM) [START_REF] Perrin | Principe du meb et microanalyse, centre de microscopie électronique à balayage et microanalyse[END_REF] (see Figure 4.1) is an electron microscopy technique capable of producing high resolution images of the surface of a sample. This kind of microscope uses an extremely fine electron beam, which is produced by thermoelectronic effect from a tungsten filament; this beam is focused at a point of the sample using an electromagnetic field. This field is generated by the coil of a condenser lens. It can, not only focus the electron beam at a single point called "cross-over", but it can also move it so that it runs through the whole sample. It is essential that the microscope column, where the electron beam spreads, remains under vacuum to avoid any deviations which could distort the measurements. Once the electron beam hits the sample, an electron-matter interaction takes place. This interaction can give birth to several types of images, and the corresponding images are shown in Figure 4.2. For our part we will focus primarily on three types of interactions:

• Backscattered electrons;

• Secondary electrons;

• X-ray.

Backscattered electrons

When the electron beam is sent in the direction of the sample, the electrons being negatively charged, they interact with the positively charged nucleus of the sample. Some of the electrons are re-emitted in a direction close to the incident direction. The sensor receiving electrons will recover high-energy electrons up to 30 keV. The amount of re-emitted electrons depends directly on the atomic number of the atoms constituting the sample. To capture the backscattered electrons, different types of sensors can be used. For example (annular semiconductors diode see Figure 4.3) axial devices are often used. In the case of semiconductors, they are divided into four quadrants A, B, C and D.

Two types of signals are conventionally used. First, the sum of the signals (A + B + C + D), which represents an average of the emission in a high solid angle. Under these conditions, we can interpret the image as a uniform light around the direction of the incident beam and as an almost purely chemical contrast. On a polished surface, the sensitivity is sufficient to detect small variations. The other kind of signal conventionally used is the signal difference between quadrants (AC, BD), the small differences in function of the backscatter emission angle are amplified. This mode enhances the topographic contrast. This is why this kind of images are used to perform homogeneity analysis. The resolution that this modality of image can reach is from one micrometer, to one tenth of a micrometer.

Secondary electrons

When part of the electrons of the electron beam reaches at the sample, a collision between the electron beam and the electrons of the sample happens. It follows that some of the electrons from the beam may transfer some of their energy to "external" electrons from the sample. By "external" we mean the farther electrons from the nucleus. This supply of energy to the "external" electrons causes their ejections of the atom. These electrons are called secondary electrons. Generally these electrons have low energy (about 50 eV) when they are ejected. Therefore, they cannot go very far and remain stuck into the surface of the sample analysed. The secondary electrons are collected with a positive electric field on a scintillator. The secondary electron efficiency received is dependent not only on the atomic number of the element observed, but also and especially on the angle between the incident beam and the sample observed. So in this type of image, the more enlightened areas correspond to areas of high efficiency, and the poor enlightened areas correspond to poor efficiency zones, one can speak of shadow area.

X-ray

The analysed sample is bombarded with a high energy electron beam about 10-40 keV. These electrons have a high energy, then they can penetrate the electron shell of atoms in the sample and excite an electron of the inner shells. The energy transfer from the electron beam to the electron in the inner shell may provoke its ejection from the atom. The ejected electron forms an electronic hole which is filled with electrons from higher energy levels. Indeed electrons from higher energy shell will gradually level down to electronic stabilization of the atom. The excitation energy of the electronic structure occurs with emission of X-ray photon, the number of X-ray photons and energy can be measured by an Energy-Dispersive Spectrometer (EDS). It is possible to characterize the different energy shells with the energy of X-ray photons and then, it is possible to characterize the electronic structure of the studied atoms. Thus we have access to information regarding the physical composition of the image. However there are two main consequences:

• Spatial resolution of the analysis is about a micron;

• The X-ray photons emitted come from a "pear" of radius and depth at about one micron (see Figure 4.4). This may pose particular problems for the analysis of small particles or complex mixtures.

Noise on SEM images

As we explained before, SEM images result from X-ray photons or electrons that are detected on various captors. Let us consider now that the level of gray of a pixel is dependent on the number of photons or electrons that the sensor receives in time interval τ , and let us write dt a smaller time interval, such as dt = τ /K, where K ∈ N. In a way we are quantifying the time interval τ . Let us consider < N (t, t+dt) > the mean number of particles ( = photons or electrons) that a sensor receives between [t, t+dt], φ the flux of particles received such that < N (t, t+dt) >= φτ /K. We will make a first hypothesis that dt is sufficiently small so that in a time interval [t, t + dt] there are just two cases:

• the sensor receives zero particles,

• the sensor receives one particle, CHAPTER 4. MULTIMODAL SCANNING ELECTRON MICROSCOPY IMAGES 107

We will write P (0) = 1q and P (1) = q respectively the probability that the sensor in a time interval [t, t + dt] receives zero particles, and one particle. Then, as one may notice, it is a binomial law of parameter q, so < N (t, t + dt) >= q. If we write the characteristic function of this law, we have:

Ψ dt (ν) = 1 -q(1 -e iν ). (4.1)
Now we have the characteristic function for a time interval of size dt, but we would like to have the characteristic function over τ . To obtain it, a second hypothesis is needed. We consider that the K intervals are independent, such that :

Ψ τ = (Ψ dt (ν)) K = 1 -q(1 -e iν ) K . (4.2)
Then we consider the asymptotic case so K → ∞, such that:

Ψ τ = lim K→∞ 1 -φτ /K(1 -e iν ) K , (4.3)
and finally we have:

Ψ τ = e (-φτ(1-e iν )) (4.4)
Therefore it is possible to recognize the characteristic function of a Poisson law of parameter φτ . Thus the number of particles in a interval τ follows a Poisson distribution of parameter φτ that can be written:

P (n) = e -φτ (φτ ) n /n!.
So contrary to hyperspectral images, where a classical assumption is that these data are corrupted by an additive Gaussian noise, here EDS data are corrupted by a Poisson noise. The processing on hyperspectral images may not necessarily work with EDS images. That is why we advice to do a variance stabilization techniques. These are techniques aiming at applying a transform Φ over the entire noisy image f corrupted by a Poisson noise so that the distribution of each pixel of the converted image, noted Φ(f ), is approximately Gaussian.

A variance stabilization technique conventionally used is the Anscombe transform [START_REF] Anscombe | The transformation of poisson, binomial and negativebinomial data[END_REF]. The transform applied to the image is:

Φ(f (x)) = 2 (f (x) + 3 8 ) (4.5)
Not only Anscomb transform stabilizes the variance of the data to 1 whatever the average value was, but also this transform will turn the data distribution to Gaussian one. After stabilizing the variance, we can perform our conventional processing, and finally we return to the data in the original space by a reverse Anscombe transform

Φ -1 (y) = y 2 2 - 3 8 (4.6) 5.1. INTRODUCTION
dicates local mean atomic number, whereas X-ray maps contains local elemental composition. Coupling a SEM with an energy dispersive spectrometer (EDS) leads to so called energy-dispersive X-ray spectrometry in the SEM (SEM-EDX). If the full X-ray spectrum is recorded on each scanned pixel, SEM-EDX produces spectral images. Typical SEM-EDX spectra are few millions pixels times few thousands of energy channels. The set of the images produced by the different SEM detectors can be seen as a multimodal image. Then, the various images can be processed independently or in a combined way. It seems also natural to consider that an improved processing would be obtained by combining the information present in the different modalities; obviously, that is true in the case where the modalities are "compatible". This process of image combination can be seen as a pratical case of the theory of information fusion. More generally, the fusion of information can be seen either as the search for optimal representation including all relevant information sources [START_REF] Piella | A general framework for multiresolution image fusion: from pixels to regions[END_REF] or as the search for algorithms making use of information from different modalities [START_REF] Wald | Some terms of reference in data fusion[END_REF].

In the context of multimodal SEM, we decided to focus on a specific problem. We aim at merging the backscattered electron images that usually have a high spatial resolution, a good signal to noise ratio but do not provide enough discriminative information to physically classify the nature of the sample, with energy-dispersive X-ray spectroscopy images that have discriminative information but a lower spatial resolution and signal to noise ratio.

This problem is similar in some ways to the so-called pansharpening [START_REF] Laben | Process for enhancing the spatial resolution of multispectral imagery using pan-sharpening[END_REF], which is well known in colour, multispectral and hyperspectral imaging. However there are also some significant differences. The first one is that in the case of classical image pansharpening, the panchromatic image (i.e., the image at the nominal spatial resolution) has a good correlation with the colour or multi/hyperspectral image. Even better, in some cases the "panchromatic image" is contained partly at some wavelengths (or linear combination of them) from the multi/hyperspectral one. In our case, the information between the two SEM modalities that are considered is not basically correlated. Moreover, we work with abundance maps extracted from the energy-dispersive spectroscopy images, since the raw EDS spectral image are of very high dimension and of very sparse nature.

As usual in multimodal SEM, we work with images perfectly registered. Although these images are spatially registered, it can be noticed that there are a number of potential artefacts which can appear when merging their underlying information sources. The origin of these artefacts is the fact that information can be structurally different, in the sense of the values of intensities, local contrast, presence of contours, etc. can be different between backscattered electron images and energy-dispersive spectroscopy. All these phenomena may degrade the quality of the fusion. A detailed list of the potential artefacts classicaly considered in image pansharpening can be found in [START_REF] Thomas | Synthesis of multispectral images to high spatial resolution: A critical review of fusion methods based on remote sensing physics[END_REF]. It seems essential to deal with these problems in order to produce image fusion without major artifacts.

Up to the best of our knowledge, the literature on SEM image fusion is inexistent. Thus to initiate our work, we decided to take inspiration from the work of fusion of information on multi/hyperspectral imaging which is much more abundant and especially on pansharpening techniques [START_REF] Vivone | A critical comparison among pansharpening algorithms[END_REF][START_REF] Loncan | Hyperspectral pansharpening: a review[END_REF]. The rest of the paper is organized as follows. After providing a description of the SEM dataset used and a summary of the most frequently used pansharpening algorithms, we introduce a new CHAPTER 5. ENHANCED EDX IMAGES BY FUSION OF MULTIMODAL SEM IMAGES USING PANSHARPENING TECHNIQUES 111 approach grounded on the bilateral filtering framework, which has been conceived in the particular case of SEM images fusion. An extensive quantitative assessment of the different algorithms is achieved to motivate discussion and conclusions.

Materials and methods

Multimodal SEM imaging.

SEM is able to produce high resolution images from the surface of a sample by means of an extremely small electron beam, which is focused at a point of the sample using the electromagnetic field of an objective lens [START_REF] Reimer | Scanning Electron Microscopy: Physics of Image Formation and Microanalysis[END_REF]. When electrons of the beam hit the specimen surface, electron-matter interactions may produce secondary particles that are detected by adequate sensors. These secondary particles can produce several types of images, as shown in Figure 1. In this work we focus primarily on two types of images :

• Images formed by backscattered electrons (BSE);

• Images formed by electron-induced X-ray detected by EDS Hence, we do not consider the Secondary electrons image, which is quite classic in SEM, due to the fact that it is less correlated with the two others modalities BSE and EDX. The different modalities acquired by the SEM provide different physical information about the sample, and consequently it can be interesting to merge them. This principle corresponds to the idea of SEM image information fusion.

SEM datasets.

SEM fusion methods discussed in this paper has been assessed using two real datasets and a simulated one.

The first dataset is composed of an EDX image of size 1024 × 768 pixels and 2024 levels of energy. Thus the image is a data cube 1024 × 768 × 2024 pixels. As an EDX image, it is naturally corrupted by a Poisson process and of extremely sparse nature, i.e., many energy levels are zero. This image is from a sample composed of iron, copper, aluminium and oxygen. This image has been acquired with a Zeiss Supra 40 SEM fitted with a 10mm 2 Brucker X-Flash 4010 EDS spectrometer. Beam energy was 15keV, probe current 0.75nA and dwell time 912µs resulting in a 717s long acquisition time. To estimate the abundance map of each of these physical elements, a Gaussian model near each peak is fitted as EDS peak shapes are very close to a Gaussian shape due to electron-hole formation statistics in the detector [START_REF] Scholze | Modelling the response function of energy dispersive x-ray spectrometers with silicon detectors[END_REF][START_REF] Bauer | Computer analysis of x-ray spectra (eds) from thin biological specimens[END_REF]. Since there are 6 detectable peaks at 6 different energies, 6 abundance maps are obtained. In the following of the paper, the multivariate image of physical abundances is called the "multispectral image". In addition to that, the BSE image of the same sample at the same spatial resolution is also available. To downscale the abundance maps, the experimental protocol that we followed is similar to the one proposed by the authors of [START_REF] Loncan | Hyperspectral pansharpening: a review[END_REF]. We have degraded the multispectral image by applying a Gaussian blur, then we have downsampled this image by a factor s = 5. The second dataset is composed of an EDX image of size 1024 × 704 pixels and 2024 levels of energy and the corresponding BSE image. This image has been acquired with a FEI Nova nanoSEM fitted with a 80mm 2 Oxford Instruments X-Max EDS spectrometer. Beam energy was 10keV, probe current 0.8nA and dwell time 3.2ms per pixel resulting in a 2307s long acquisition time. This image is from a sample composed of aluminum, oxygen, vanadium and nickel. In this case, 4 abundance maps were obtained using the same model and algorithm. Downscaled images of the abundance maps image were produced too.

In order to assess the accuracy of the algorithms, EDX images have been simulated by the Monte-Carlo method by a dedicated software based on the PENE-LOPE package [START_REF] Sempau | An algorithm for monte carlo simulation of coupled electron-photon transport[END_REF][START_REF] Sempau | Experimental benchmarks of the monte carlo code penelope[END_REF][START_REF] Salvat | Penelope-2006: A code system for monte carlo simulation of electron and photon transport[END_REF]. First, the response of the EDS detector of the Nova NanoSEM has been characterized by measuring the full width at half maximum (FWHM) for different peak energies. The dependence on peak energy E of the standard deviation of the Gaussian peaks σ(E) was obtained by least-square fitting and gave :

σ(E) = √ 0.4403E + 337.04, (

with E in eV. The efficiency (E) of the detector was modeled from its geometry following the same approach proposed by Limandri et al. [START_REF] Limandri | Experimental study of the efficiency of a sdd x-ray detector by means of pixe spectra[END_REF]. Analog electrons trajectories (C1 = C2 = 0; W CC = W CR = 0) were generated with the PENE-CHAPTER 5. ENHANCED EDX IMAGES BY FUSION OF MULTIMODAL SEM IMAGES USING PANSHARPENING TECHNIQUES 113 LOPE package. As X-ray emission is an unusual process, both characteristic and Bremsstrahlung photon emission were enhanced by interaction forcing by a factor F with the help of the build-in functions of PENELOPE. Sample and detector geometries were handled by the PENGEOM package. The detector is annular at a 35 • take-off angle to the surface. Any photon hitting the detector was considered detected with a probability (E). Each detected photon of energy E was registered in an energy channel distributed following a Gaussian law of average E and standard deviation σ(E).

State-of-the-art

Before presenting the different image fusion methods, let us introduce the notation used in the following, which is inspired from [START_REF] Loncan | Hyperspectral pansharpening: a review[END_REF]. First, from a mathematical viewpoint, a multispectral image (abundance maps) is considered as a function HS defined by

HS : E → R D x → v i
where D is the number of abundance maps and E is the image domain (support space of pixels). This multivariate image can be also seen as set of D grey-scale images.

We note HS a multispectral image at a low resolution. Let HS ∈ R N 1 ×N 2 ×D be a interpolated multispectral image whose spatial dimensions are N 1 , N 2 , which in our case corresponds to BSE image dimensions. Let R ∈ R N 1 ×N 2 be just the BSE image. We denote by HS ∈ R N 1 ×N 2 ×D the multispectral image of the EDS abundances enhanced with the BSE image information, where HS k is the k-th abundance map of image HS.

There are essentially three families of pansharpening techniques which are detailed as follows.

Component substitution methods (CS)

The purpose of information fusion techniques is to find the function φ satisfying

HS = φ( HS, R),
where HS is "optimal" in a certain sense. The particular notion of optimality is precised below. Component substitution methods (CS) are based on the projection of the image HS into another space, separating the spatial information from the spectral one. The spatial information, i.e., contrast and contours between different objects, is often concentrated in a single grey-scale image. The main step of CS consists in replacing the image containing the spatial information by the spatially high resoluted image R. It is based on the assumption that the image containing only spatial information is highly correlated with R. Once the corrections are done, the data are projected back to the initial space. This approach is global and therefore corrects all pixels with a single rule. The strength of this technique is its speed, however it depends on the fact that the image containing the spatial information of HS which must be similar in range and intensity distribution to R. If not, the merged image may have strong distortions. In our particular case, the BSE image has information which is not comparable with EDS images, that explains the artefacts which appears in our results, see Section 5.

The problem can be formulated mathematically as follows [START_REF] Loncan | Hyperspectral pansharpening: a review[END_REF]:

HS k = HS k + g k (R -I L ), (5.2)
where g k corresponds to a corrective coefficient for the band k and I L is the image containing all high resolution spatial information, i.e., typically

I L = D i=1 w i HS i , (5.3)
where w i is a coefficient depending on the degree of the spatial information of the spectral band i. Now let us present in a more detail way the most popular CS techniques.

CS by PCA.

Principal Component Analysis (PCA) [START_REF] Pearson | Liii. on lines and planes of closest fit to systems of points in space[END_REF][START_REF] Hotelling | Analysis of a complex of statistical variables into principal components[END_REF][START_REF] Chavez | Comparison of three different methods to merge multiresolution and multispectral data-landsat tm and spot panchromatic[END_REF] is a well-known dimensionality reduction technique, which aims at finding an orthonormal basis that maximizes the variance of the data. Thus, the data projected into this space summarize the statistically significant information. The fundamental assumption of this pansharpening technique is that the first component focuses all relevant spatial information while the other components contain secondary spatial and spectral contrast features. Using the CS model, the first principal component of HS is I L , and w i is the first column of the inverse transform. Note that this is not exactly the R which is used 5.3. STATE-OF-THE-ART in the Eq. ( 5.2), since the version of R should be histogram equalized with respect to I L .

CS by Gram-Schmidt decomposition (GS).

The Gram-Schmidt (GS) technique is based on an orthogonal decomposition, which was invented by Kodak [START_REF] Laben | Process for enhancing the spatial resolution of multispectral imagery using pan-sharpening[END_REF]. First, instead of considering only the cube HS, one starts with the tensor composed of the concatenation of HS and R, thus of D + 1 bands. Then, an orthogonal decomposition is performed on this new tensor. We find the vector of the basis that corresponds to the spatial information and it is replaced by R. Finally a reverse procedure is done to return to the initial representation space. The components of the Gram-Schmidt basis being orthogonal, the spatial information is expected to be orthogonal to all the other information sources. That is why the spatial information is, in theory, gathered in a single vector. In practice this is not always the case, especially if R contains dense information. Moreover this technique can be problematic in the case of multimodal SEM images since a significant orthogonal basis from EDX + BSE cannot be easily obtained without adding a kind of sparsity constraint.

Multiresolution analysis methods (MRA)

MRA techniques are founded on the application of a low pass filter to R, sometimes under different resolutions. Then, details are injected on HS thanks to the residue of R L and R, which represents the high frequences of R. A mathematical formulation of this type of technique is written as [START_REF] Ranchin | Fusion of high spatial and spectral resolution images: the arsis concept and its implementation[END_REF][START_REF] Tu | A new look at ihs-like image fusion methods[END_REF]:

HS k = HS k + G k ⊗ (R -R L ), (5.4)
where G k corresponds to a corrective coefficient matrix for the band k and ⊗ is the multiplication term by term. Moreover, as mentioned above, R -R L represents the details that are injected in the low resolution image. Such family of techniques depends mainly on the type of decomposition performed for R L and the diagonal matrix of gains G k . Different kind of filters can be used for R L : Gaussian filters, wavelets, mathematical morphology operators, etc. There are therefore many possibilities. Two techniques are now detailed.

MRA by Smoothing Filter-based Intensity Modulation (SFIM).

The technique [START_REF] Liu | Smoothing filter-based intensity modulation: a spectral preserve image fusion technique for improving spatial details[END_REF] involves the use of a single low pass filter, noted H lp , applied to R to get R L . The enhancement is then obtained as:

HS k = HS k + G k ⊗ (R -R * H lp ), (5.5)
with

G k = HS k R * H lp , (5.6)
where is the division term by term. To achieve the low pass filter, a simple average filter is typically used. The formula can be simplified thanks to the coefficient matrix to obtain:

HS k = R ⊗ ( HS k R l ). (5.7)
CHAPTER 5. ENHANCED EDX IMAGES BY FUSION OF MULTIMODAL SEM IMAGES USING PANSHARPENING TECHNIQUES 117 Thus it involves that the information of R was modulated by the ratio of the image HS k and R l . This allows integrating the contrast present in R without creating missing objects in the multispectral image HS. Due to this pixel by pixel multiplication, there can be an issue on the low values of R l , which requires to manage the dynamic of the image, such that it does not exceed a given range.

MRA by Laplacian Pyramid (MTF-GLP).

This technique [START_REF] Burt | The laplacian pyramid as a compact image code[END_REF][START_REF] Mallat | A theory for multiresolution signal decomposition: the wavelet representation[END_REF] has some similarities with the previous one. However, this time a multiple low-pass filter H reso i at various resolutions (i.e., at various scales) is used. More precisely, the low-pass filters are typically Gaussian convolutions. Missing information at a given scale i is injected into HS reso i+1 thanks to information from R at the resolution i. In our framework, the BSE image R and the multispectral abundance maps image of the EDS image are both in a multi-resolution structure pyramid, as illustrated in Figure 5.5. Details are injected at each resolution, until the nominal resolution of the BSE image is reached.

An hybrid method: Guided PCA

This technique [START_REF] Liao | Processing of multiresolution thermal hyperspectral and digital color data: Outcome of the 2014 ieee grss data fusion contest[END_REF][START_REF] Kang | Spectral-spatial hyperspectral image classification with edge-preserving filtering[END_REF] consists in first doing a PCA on the multispectral image HS. Then, on the d first principal components, a guided filter [START_REF] He | Guided image filtering[END_REF] is applied and nothing is done on the remaining Dd bands. The rationale is based on the fact that the remaining bands correspond to "noise", and it would be useless to enhance the resolution of noise. Then a reversed PCA transform is made to obtain the original multispectral representation. The high resolution information from R is included by the guided filter. For more details about this technique, see [START_REF] Liao | Processing of multiresolution thermal hyperspectral and digital color data: Outcome of the 2014 ieee grss data fusion contest[END_REF][START_REF] Kang | Spectral-spatial hyperspectral image classification with edge-preserving filtering[END_REF].

Fusion of SEM information by Abundance Guided

Bilateral Filter (AGB)

As we have discussed above, when a pansharpening technique is used in the image fusion context, the starting point is a low spatial resolution image HS, together with a high resolution "panchromatic" image R. Usually, the first step consists in upsampling HS, thanks to a basic bi-cubic interpolation to obtain the image HS. Then, on this image at the nominal scale, the image enhancement is done, where the information is corrected by means of a particular pansharpening algorithm. The main issues with the previous techniques are that the CS approaches do not take into account the differences between the various components of the multimodal image, that is, the differences between the abundance maps themselves and the differences with the panchromatic image. In our case, there is a few correlation between the different EDX abundance maps. Moreover the MRA approaches do not consider the spectral link between the different abundance maps, since there is a finite quantity of material for each pixel (i.e., the sum of abundances at a given pixel is equal to 1).

In order to address these drawbacks, we propose an interpolation technique called abundance guided bilateral filter (AGB) by considering the relationships between the abundance maps. Thus, weights used in the interpolation would depend on both R and HS. Our approach of interpolation uses bilateral filter [START_REF] Tomasi | Bilateral filtering for gray and color images[END_REF] and more exactly its cross version [START_REF] Eisemann | Flash photography enhancement via intrinsic relighting[END_REF]. By the way, the bilateral filter has already been used on other super-resolution algorithms. Inspired by these works, we have chosen to improve the interpolation process which is a scaling process.

A bilateral filter is a nonlinear, edge-preserving denoising/regularizing operator for images. The intensity value at each pixel in an image is replaced by a weighted average of intensity values from nearby pixels. This weight can be based on a Gaussian distribution. Crucially, the weights depend not only on Euclidean distance between pixels on the grid, but also on their intensity (or more general radiometric) differences. Thanks to the fact that it uses spatial and range information, it preserves the edges, this is the reason why bilateral filter is used in super-resolution.

Formally, the joint bilateral filter of an image I guided by an image F is defined as :

I * (x) = 1 W p (x) x i ∈E I(x i )k(x, x i ), (5.8)
where the kernel weights are given by

k(x, x i ) = f r ( F (x i ) -F (x) )g s ( x i -x ),
and where the normalization term is just given as:

W p (x) = x i ∈E k(x, x i ),
such that f r and g s represent respectively the range (or spectral) and spatial kernels of parameters r and s. To simplify, we have chosen a Gaussian function for both kernels. However in our case we would like to consider the link between the different abundance maps. To handle this relationship, we need to define a guided function as a global criterion. Hence, we used the level of mixing of pixels, which will be represented by an order map: an image of ordered levels of intensity. In a way the order map has a low value if the pixel contains a mixture of various different materials and a high value if it is almost pure, so it contains contribution of few materials.

CHAPTER 5. ENHANCED EDX IMAGES BY FUSION OF MULTIMODAL SEM IMAGES USING PANSHARPENING TECHNIQUES 119 Let us precise how this order map is computed using the possible alternatives. The EDS spectrum at a position i of the image is a vector v i which can be written as:

v i = D r=1 a r,i m r + n i , (5.9)
where {m r } D r=1 represent the set of D endmembers (spectral signature of the material), a r,i the abundance at vector i of each endmember r, and n i an additive noise. This last term can be neglected. The nonnegative coefficients a r we consider are convex combination such that D r=1 a r = 1. By using this additional constraint, it is guaranteed to work on a (D -1)-simplex.

We have proposed in [START_REF] Franchi | Ordering on the probability simplex of endmembers for hyperspectral morphological image processing[END_REF] different techniques to calculate such order map on abundance map images. We adopt here an approach based on the notion of majorization [START_REF] Olkin | Inequalities: theory of majorization and its applications[END_REF], which is a technique for ordering vectors of same sum. Let us consider two vectors c = (c 1 , . . . , c n ) ∈ R n and d = (d 1 , . . . , d n ) ∈ R n , then we say that C weakly majorizes D, written c w d, if and only if

k i=1 c ↓ i ≥ k i=1 d ↓ i , ∀k ∈ [1, n] n i=1 c i = n i=1 d i (5.10)
where c ↓ i and d ↓ i represent respectively the coordinates of C and D sorted in descending order. Majorization is not a partial order, since c d and d c do not imply c = d, it only implies that the components of each vector are equal, but not necessarily in the same order.

Let us define a majorization-like partial order adapted to the abundances. A permutation τ i of the coordinates of the vectors v i in the simplex is applied such that they are sorted in descending order. The majorization-like order ≤ maj is defined as

v i ≤ maj v j ⇔                  a τ -1 i (1),i < a τ -1 j (
1),j or a τ -1 i (1),i = a τ -1 j (1),j and a τ -1 i (2),i < a τ -1 j (2),j or . . .

a τ -1 i (1),i = a τ -1 j (
1),j and . . . and

a τ -1 i (R),i ≤ a τ -1 j (R),j
This order map between pixels brings a global information: the materials entropy. Since a pixel that has a high value with this order is less mixed than the other one. Let us write O the order map of the abundance images from the EDS image. Figure 5.6 provides the corresponding O image for dataset 1 at the resolution of HS, together the BSE image R at the nominal resolution. Let us write O the order map image at the nominal resolution, where the value of all the missing pixels is put to zero.

We can now introduce the expression of the SEM fusion using the AGB filter: where the AGB kernel is written as the product of two kernels

HS(x) = 1 W p (x) x i ∈E HS(x i ) k(x, x i ) ( 5
k(x, x i ) = k space (x, x i )k EDX+BSE (x, x i ) (5.12)
given by

k space (x, x i ) = g s ( x i -x )M ( O(x i )), k EDX+BSE (x, x i ) = g s (| O(x i ) -O(x)|)g s (|R(x i ) -R(x)|),
with a typical normalization term:

W p (x) = x i ∈E k(x, x i ), (5.13)
M is a mask that allows us to consider just the data on the values of abundance images available at the low resolution:

M ( O(x i )) = 1, if O(x i ) = 0 0, otherwise (5.14)
and finally, g s is a Gaussian function of scale parameter s.

Fusion of SEM information by Bilateral Guided Morphological Filter (BGM)

The goal now is to use mathematical morphology operators to build a nonlinear interpolation. This time, instead of using a classical average on the kernel, to approximate the missing data, we will use an alternative mean called the L ∞ barycenter, which can be seen related to mathematical morphology operators. Let us first consider a grey-scale image f : E → T, T ⊂ R. The two basic operations in morphology are the grey-level erosion and the grey-level dilation whose definition respectively are:

ε b (f )(x) = inf h∈E (f (x -h) -b(h)),
(5.15)

δ b (f )(x) = sup h∈E (f (x -h) + b(h)), (5.16)
where b is a structuring function, which introduces the effect of the operators by the geometry of its support as well as its weights. Here we consider a structuring function with weights related to the bilateral kernel.

In order to compute the dilation and erosion, we obviously need a partial order for the supremum and infimum. Thus, we can use the order previously develop, represented on the image space by the order map O. Using this order, we can naturally define an erosion and dilation on the multivariate image of abundances as follows:

ε b (HS)(x) = HS(x a ), (5.17) where

x a = arg min h∈E (O(x -h) -b(h)), (5.18
)

δ b (HS)(x) = HS(x b ), (5.19)
where

x b = arg max h∈E (O(x -h) + b(h)). (5.20)
Moreover, in our case, the dilation and erosion that we use involve the concept of changing the scale of the images. This notion is related is linked with the theory of morphological wavelets. We are not going to recall all the concepts that can be found in [START_REF] Heijmans | Nonlinear multiresolution signal decomposition schemes. ii. morphological wavelets[END_REF][START_REF] Goutsias | Nonlinear multiresolution signal decomposition schemes. i. morphological pyramids[END_REF], but at least to basic notions required for the rest.

In signal processing, the Shannon theorem states that the sampling of a signal, that is a discrete representation by a set of regularly sampled values of the signal, requires a frequency greater than Nyquist frequencies. In the case that this theorem is not respected the obtained signal may suffer from aliasing effects. This is the reason why Haralick et al. [START_REF] Haralick | The digital morphological sampling theorem[END_REF] defined a sampling condition for image processing that we are going to introduce.

Theorem 22 Let S ∈ Z 2 be a subset representing the sampling set of the grey-scale image f defined on the domain E ⊂ Z 2 , and let K ∈ Z 2 be a flat structuring element (SE), i.e., a set defined at the origin. This structuring element is large enough to verify the following sampling condition if:

δ K (S)(x) = E.
(5.21)

This condition is important to reconstruct the signal since it means that the structuring element is large enough to cover all E. According to [START_REF] Heijmans | Morphological sampling[END_REF] the sampling condition implies that the sampling distance must be less than half the distance of the structuring element K, that can be rewritten mathematically as:

x ∈ K y =⇒ K x ∩ K y ∩ S = ∅ K ∩ S = {0, 0}
(5.22) MORPHOLOGICAL FILTER (BGM)

where K x = {k + x|k ∈ K}. Based on this assumption, Heijmans and Toet defined another general sampling strategy that we are going to review. Let E and S be two subsets Z 2 , E representing the support space of f and S the one of the sampled version of f . Let K : S -→ P(E) be a mapping where P(E) represents the power space of all the subsets of E, where ∀s ∈ S, K(s) = {k + s|k ∈ K} that means that it is the translate of K along s. Then it is also possible to define the dual mapping as:

K * : E -→ P(S) K * = {s ∈ S|x ∈ K(s)} (5.23)
Thanks to this definition, it follows that:

K * * = K ∪ s∈S K(s) = E if and only if ∀x ∈ E, K * (x) = ∅.
That means that K covers E. Let us consider f : E → F a grey-scale image and f S : S → F S its sampling, and let us consider fun(F ) the set of grey-scale images and fun(F S ) the set of grey-scale sampled images. It is now possible to define the dilation sampling operator:

σ K (f ) : fun(F ) -→ fun(F S ) σ K (f )(s) = sup{f (x)|x ∈ K(s)} (5.24)
and the erosion sampling operator: We can prove that σ and σ are adjoint operators using the standard procedure

ν K (f ) : fun(F ) -→ fun(F S ) ν K (f )(s) = inf{f (x)|x ∈ K(s)} (
σ K (f ) ≤ g S ⇔ ∀s ∈ S, sup{f (x)|x ∈ K(s)} ≤ g S (s) ⇔ ∀s ∈ S, ∀x ∈ K(s), f (x) ≤ g S (s) ⇔ ∀x ∈ E, ∀s ∈ K * (x), f (x) ≤ g S (s) ⇔ ∀x ∈ E, f (x) ≤ inf{g S (s)|s ∈ K * (x)} ⇔ f ≤ σK (g S )
Identically we have that ν and ν are adjoint too. Therefore, according to the fundamental theorem of mathematical morphology [START_REF] Heijmans | Mathematical morphology and its applications to image and signal processing[END_REF] , it is possible to build new operators by composition of the adjoint operators: where ρ is a closing and η is an opening. Heijmans and Toet named σ and ν the sampling operators while ρ and η are called the reconstruction operators. We follow their terminology. These sampling/reconstruction operators are at the basis of the morphological wavelet theory. We notice that these operators used a fixed structuring element K depending on just one image. By doing it, they do not use the intrinsic property of the images and so the sampling is not image-adaptive, moreover they cannot use the information from a second image. One would like that the sampling depends on both the original image and another image to have a morphological wavelet operators useful for image fusion.

ρ K = σK σ K (5.28) η K = ν K ν K (5.
Let us write Õa the order map image at the nominal resolution (i.e., the resolution of image R), where the value to all the missing pixels is put to ⊥. The order map image Õb is that where the missing values are put to , where ⊥ is smaller than all the values of O and is higher than all the value of O. We introduce now two new operators: σK (HS)(x) = HS(x b ), (5.30) where

x a = arg min h∈E ( Õa (x -h) -K x (h)), νK (HS)(x) = HS(x b ), (5.31)
where

x b = arg max h∈E ( Õb (x -h) + K x (h)),
This definition implies that σK (HS) and νK (HS) are defined at the resolution of R.

In order to improve the accuracy of our morphological operators we choose to use bilateral structuring function, inspired by the theory developed in [START_REF] Landström | An approach to adaptive quadratic structuring functions based on the local structure tensor[END_REF][START_REF] Angulo | Morphological bilateral filtering[END_REF]. The bilateral structuring function depends on x and its neighbourhood. However since we want to add the information from R, it should be also depending on the latter:

K x (h) = 1 Wp g s ( Õa (h) -Õa (x) )g s ( h -x )g s ( R(h) -R(x) )M ( Õ(h)) if h ∈ Ω x -∞
elsewhere (5.32) where the normalization term is:

W p = max h∈Ωx g s ( Õa (h) -Õa (x) )g s ( h -x )g s ( R(h) -R(x) )M ( Õ(h)) . (5.33)
We wrote the structuring function for the case of the locally adaptive erosion; in order to compute the corresponding locally adaptive dilation, one just has to replace Õa by Õb on the two previous formulae.

By means of these operators, it is possible to increase the scale of an image based on its own information but also the information from R. Using only one of the operators, the fused image may have a distortion by skewing the values to the high or low values. In order to correct this effect, we use as an enhanced image the mean between the adaptive upsampling erosion and dilation, i.e., ĤS = σK (HS) + νK (HS) 2 .

(5.34)

This operator is just a kind of locally adaptive L ∞ barycenter.

Evaluation criteria of pansharpening algorithms

The criteria used to evaluate the quality of the merged information are those conventionally considered in the context of the pansharpening literature. More precisely, we dealt with the following criteria:

C1: The spectral distortion between the enhanced multispectral image and the real multispectral image at the nominal resolution should be as small as possible. Or, in other terms, we would like to find the same materials for each pixel as the original image at high resolution;

C2: The spatial distortion between the enhanced multispectral image and the real one should not be too high.

Many alternative metrics can be used to quantify these two criteria [START_REF] Loncan | Hyperspectral pansharpening: a review[END_REF]. Let us formaly precise those that we have used. We write HS the multispectral abundance EDS image at the same resolution that R, that has been provided by the sensor at high resolution. In a way it is the ground truth such that our enhanced images HS should be compared to HS. Moreover, we denote by M ∈ M n,D (R) and M ∈ M n,D (R) the two matrices representing respectively HS and HS, where n is the total number of pixels (i.e., n = N 1 × N 2 ), and D the number of abundance maps.

We write by M i,: , ∀i ∈ [1, n], a spectra of HS, and M :,k ∀k ∈ [1, D] a map of HS.

We have now all the notation to introduce the four parameters. Cross correlation (CC) is a measure that evaluate the spatial distortion defined as

CC( HS, HS) = 1 D D k=1
CCS( M :,k , M :,k ), (5.35) where

CCS( M :,k , M :,k ) = ( n i=1 M i,k -µ M :,k )( n i=1 M i,k -µ M :,k ) n i=1 ( M i,k -µ M :,k ) 2 n i=1 ( M i,k -µ M :,k ) 2
, with µ M :,k = n -1 n i=1 M i,k being the empirical mean. The CC is optimal when it is close to 1. where < •, • > is the Euclidean inner product of vectors associated to the norm L 2 , and where • is the norm L 2 of vectors. The SAM is optimal when it is near to 0.

Spectral Angle

Root mean squared error (RMSE) measures the mean residual error of fusion and is defined as:

RM SE( HS, HS) = M -M F n.D , (5.37)
where • F is the Frobenius norm of a matrix A, i.e., A F = trace(AA t ). The RMSE is optimal when it is near to zero.

Synthetic adimensional global error (ERGAS) offers an overall measure of the quality of an enhanced image. It is given by the expression:

ERGAS( HS, HS) = 100d D k=1   RM SE( M :,k , M :,k ) µ M :,k   2 ,
where d is the ratio between the linear resolution of the BSE image R and the abundances EDS image HS, i.e., d = R-linear spatial resolution HS-linear spatial resolution .

The ERGAS is optimal when close to 0.

We note that C1 and C2 are not enough to assess in our case the quality of the enhancement. Indeed the BSE images give not just access to the high resolution spatial information of the physical objects, but it provides other kind of information that can in theory allow us to have better result than just the EDX image at high resolution, such as the topographic shape of the sample. Thus we need to introduce a new criterion measuring the amount of information which is injected in the merged image: C3: The injected information. Since the information provided by the various modalities can be relatively different, we would like to inject the useful one to improve segmentation and characterisation of the EDX image.

In order to quantify this criterion, we proposed to use the cross correlation between the norm of gradient of the images. To calculate such norm of the image gradient, different techniques can be used. In our studies, we computed the morphological gradient [START_REF] Soille | Morphological image analysis: Principles and applications[END_REF]. More precisely, we have:

Cross correlation gradient (CC g ) is a measure that evaluate the spatial distortion defined as

CC g (R, HS) = 1 D D k=1 CCS(Rg, Mg :,k ),
where Rg represents the image gradient of the BSE image converted into a vector and Mg :,k is the image gradient of the enhanced abundance map at the nominal scale converted also into a vector.

An optimal enhancement method would be a compromise between the criterion C3 and the C1 and C2.

Evaluation on dataset 1

The EDS abundance maps of dataset 1 are provided in Figure 5.7. The results of the enhanced abundances HS, obtained by the different techniques are provided in Figures 5.8,5.9,5.10,and 5.11 where the abundances are visualized as RGB color images. Quantitative results for this case according to the five measures are presented in Table 5.1.

From these results, we note that CS techniques are good to inject R on HS, while MRA techniques are good to increase the resolution. In addition, AGB, GFPCA and MTF-GLP present both a good compromise between these criteria. The BGM technique seems to be good to inject the information expect in the pixels with a high gradient.

Evaluation on dataset 2

In the case of this multimodal SEM dataset, which has a higher level of noise, we try to inject the BSE image R, depicted in Figure 5.12, to increase the resolution of the EDS abundance maps provided in Figure 5.13, and also to denoise the images. The results for our ABG method are given in Figure 5.14. For comparaison, Figures 5.15 and 5.16 provide a visualization of the results obtained using the other methods. Quantitative results are given in Table 5.2, which lead to similar conclusions as for dataset 1.

Evaluation on simulated dataset

By means of the Monte Carlo software based on PENELOPE, we calculated various simulated datasets with different levels of noise and different resolutions. So first we have a dataset represented in Figure 5.2. It represents the sample at the full resolution without noise, thus it corresponds to the ground truth. The size of this image is 1024×1024×2048. Then we simulate other EDX images of the same sample at resolution 4 times smaller than the original. The size of these images is therefore 256×256×2048. First we extract the 7 abundance maps of these images and perform the fusion techniques. The results of the different algorithms are gathered in Table 5.3. Since we know the ground truth, for each pixel the most important materials that are present are known (see Figure 5.17). Hence after having performed the fusion techniques, for each pixel we extract the most significant material. This imply that, each pixel is represented by a number between 1 and 7. Then these images can be seen as the results of a classification. To evaluate the performance we used the overall accuracy (OA) which is defined as the percentage of pixels well classified on the whole image, and the average accuracy (AA) which represents mean of the accuracy of each class. These grey scale images are represented in Figure 5.17.

One can see the results with the different levels of noise in Table 5.4. Our proposed ABG method and GFPCA give the best accuracies, even with high noise levels.

Conclusion

In this article, we have conducted a qualitative and quantitative evaluation of different algorithms of fusion of information applied on SEM images. Thanks to these algorithms we merged EDX and backscattered electrons SEM images. We compared global and local state of the art techniques and proposed two innovative one based on a guided bilateral filter, and a morphological approach. We also wrote a program based on the PENELOPE package to simulate large SEM-EDX spectral maps. Thanks to these simulated and real SEM images, we evaluated the performances of the different techniques according to criteria defined for pansharpening problem. We also proposed one new criterion that is more adapted for the SEM segmentation problem. The accuracy of the pansharpening methods was assessed, and our proposed method was always found to be one of the best, together with the GFPCA. Our proposed algorithm performs well even on very noisy EDX images. Thanks to this fusion of information, we can increase the speed of the imaging process, since we need fewer pixels. A way to improve the compression of this technique might be to consider pixels at given location, where the location is guided thanks to the backscattered image. By doing that, we would first take the backscattered image then have a look at the keypoints localizations, and then take the EDX image just at those points. Another way to improve this technique might be to consider the secondary electrons image. In addition, thanks to the developed formalism, it can easily be taken into account. In a broader perspective, the proposed guided bilateral filter could be used for other multispectral image pansharpening than SEM based ones.

Chapter 6

Spatial Regression for Image Pansharpening: Application to Multimodal SEM Image Fusion

Abstract

In this chapter, we compare ordinary kernel regression, kriging and Gaussian processes and apply such methods to multimodal multispectral SEM image fusion. The kernel regression that we use is the technique presented in the previous chapter, called Abundance Bilateral Guided (ABG). We compare it mathematically and experimentally to ordinary kriging. In addition we propose a way to perform image fusion of information by means of the ordinary kriging and show the interest of the approach for image pansharpening. By combining the different SEM modalities, we increase the resolution of the EDX multispectral image, and provide more details on the physicochemical composition of samples obtained from the enhanced EDX.

Résumé

Le microscope électronique à balayage (MEB) permet d'acquérir des images à partir d'un échantillon donné en utilisant différentes modalités. Le but de ce chapitre est d'analyser l'intérêt de la fusion de l'information pour améliorer les images acquises par MEB. Nous avons mis en oeuvre différentes techniques de fusion de l'information des images, basées dans ce chapitre sur la théorie de la régression spatiale. Ces solutions ont été testées sur quelques jeux de données réelles.

Introduction

In this chapter, we propose an original approach of spatial fusion of information able to increase the resolution of multispectral images. This technique is called pan sharpening. Pan sharpening aim at merging different spatial and spectral information. A huge variety of techniques exists in remote sensing [START_REF] Laben | Process for enhancing the spatial resolution of multispectral imagery using pan-sharpening[END_REF], however in Scanning Electron Microscopy these kind of techniques are 6.1. INTRODUCTION not so used. We express in this chapter this problem as a problem of spatial regression. To perform this regression we consider different techniques. We compare them mathematically and propose one innovative technique based on regularised kriging. We note that a solution based on Parzen windows was introduced in the previous chapter. So we present the mathematical aspect of this solution, and compare it to our solution based on kriging. Kriging [START_REF] Matheron | The intrinsic random functions and their applications[END_REF][START_REF] Chiles | Geostatistics: modeling spatial uncertainty[END_REF][START_REF] Cressie | Statistics for spatial data[END_REF][START_REF] Matheron | Principles of geostatistics[END_REF][START_REF] Cressie | The origins of kriging[END_REF] is an interpolation technique used originally in geosciences for the evaluation of minerals deposit spatial repartition. This method was developed by Georges Matheron departing from some studies from Daniel G. Krige. The technique based on kriging is inspired by the innovative work [START_REF] Dabo-Niang | Combining kriging, multispectral and multimodal microscopy to resolve malaria-infected erythrocyte contents[END_REF] where the author first performed a principal component analysis followed by kriging applied on each principal component. By combining this information with other modalities they were able to recognize if a patient was infected by the malaria or not. We note also the innovative work proposed by [START_REF] Li | Learning local gaussian process regression for image super-resolution[END_REF], where the authors use local gaussian process regression techniques to increase the resolution of an image. Given an image the authors randomly extract patches from a blurred version of the image and on its corresponding high-resolution version. Then given pairs of patches they learn thanks to a technique based on multiple Gaussian process regressions how to transform a low resolution patch into a high resolution one. Even if this technique is quite interesting, the Gaussian process regression they used does not take into consideration the full spatial information of images.

That is why, here we propose to use a local ordinary kriging where the coefficients depend on the information of the other modalities we want to inject. Then we are able to reconstruct the missing information while injecting other modalities information. To present our work we first introduce the kernel regression, then gaussian process regression, then ordinary kriging. We provide results of these techniques on the samples presented in the previous chapter.

Notations

We consider that an image f of dimension D with support space of pixels E ⊂ Z 2 can be defined as a function:

f := E → F ⊂ R D x → f (x) (6.1)
In the sequel, we limit ourselves to two kinds of images from multimodal SEM. We first consider the backscattered electron (BSE) image, which is denoted R. It is a gray scale image, i.e., D = 1, of spatial dimensions N 1 × N 2 . We also consider a multispectral image HS which is in our case the result of the abundance maps, at a lower resolution, of dimension n 1 × n 2 × D, where D is the number of abundance maps. We additionally denote HS the abundance map multispectral image at the nominal spatial resolution: size N 1 × N 2 × D. The goal of image fusion is just to increase the spatial resolution from n 1 × n 2 to N 1 × N 2 .

We consider that an image g is the realization of a random function G. Let us consider a probability space (Ω, A, P ) and a domain D ∈ E. A random field, or random function, on the spatial domain D with values in F is a function of two variables, denoted Z(x, ω). For each x ∈ D, G(x, :) : ω → G(x, ω) is a random variable on (Ω, A, P ). Moreover, for each ω 0 ∈ Ω, G(: ω 0 ) : x → G(x, ω 0 ) is a CHAPTER 6. SPATIAL REGRESSION FOR IMAGE PANSHARPENING: APPLICATION TO MULTIMODAL SEM IMAGE FUSION 141 function of D → E. We write G(x) the random function at the position x, and g(x) the regionalized function.

We note that in our case we deal with image data hs that are abundance map. That means that at a given position x, we have D i=1 hs i (x) = 1, together with hs i (x) ≥ 0, ∀x ∈ E, ∀i. This means that the data are on a convex set of the positive orthant of R D , called D -1 simplex. From a physical viewpoint, hs i (x) represents a quantity (between 0 and 1) of material i at pixel x.

Spatial kernel regression

Let us consider that the image hs is composed of N spatial positions N 1 = {x i |i ∈ [1, N ]}, such that we have the value of hs in just n spatial positions N 2 = {x i |i ∈ [1, n]}. Our goal is to estimate the value of hs in the remaining Nn spatial positions

N 3 = {x i |i ∈ [1, N ]/[1, n]}.
Or, more formally, given a vector of inputs N 2 = (x 1 , x 2 , ..., x n ) t , we will predict the output hs for new positions x * . The estimated predictor is denoted hs(x * ).

Let us consider that the available information that we have on each observation is most of the time corrupted by Gaussian noise of zeros mean and standard deviation σ. Hence, we can write the model in the form:

hs(x i ) = f (x i ) + N (O, σ) (6.2)
where f is the true information. By considering that the different data points (x i , hs(x i )) are independent and identically distributed (i.i.d.), we can evaluate the likelihood as follows:

P(HS = hs(N 2 )|N 2 , f, σ) = n i=1 P(HS = hs(x i )|x i , f, σ), (6.3)
with hs(N 2 ) = (hs(x 1 ), . . . , hs(x n )) t .

P(HS

= hs(N 2 )|X, f, σ) = n i=1 N (hs(x i )|x i , f, σ) (6.4)
P(HS = hs(N 2 )|N 2 , f, σ)) = 1 √ 2π n σ n .e 1 2σ 2 f (N 2 )-hs(N 2 ) 2 (6.5)
We seek a function f (x) for predicting hs(x) given values of the input N 2 . However, the optimal function should not depend on N 2 . This leads us to a criterion for choosing f , in the sense that the optimal function f should minimize:

R(f ) = E hs,x ( f (x) -hs 2 ) = E x E hs|x ( f (x) -hs 2 ), (6.6)
which can be rewritten as:

f (x) = arg min E HS|X=x ( f (x) -hs 2 ), (6.7)
6.2. SPATIAL KERNEL REGRESSION and its solution is just:

f (x) = E HS|X=x (hs), (6.8)
which corresponds to the conditional expectation, also called the regression function.

Then the best prediction needs P(HS|X = x). However, we do not have access to this probability and the goal here is to show how to evaluate it thanks to kernel density estimation also known as Parzen-Rosenblatt window method [START_REF] Rosenblatt | Remarks on some nonparametric estimates of a density function[END_REF].

We consider a model of the data which is similar to the one introduced in [START_REF] Bishop | Pattern recognition[END_REF]. Let us write T = (HS, X) the joint random variable. For the sake of simplicity, we write p(t) the corresponding joint distribution with t ∈ R D+2 . Let us take a (small) region R ⊂ R D+2 containing t. The probability of this region is just:

P = R p(t)dt.
(6.9)

We suppose that out of the n points there are K points that fall into R, such that the K points follow a binomial distribution, i.e.,

Bin(K|n

, P) = n! K!(n -K)! P K (1 -P) n-K . (6.10)
If we want the mean number of points falling in R to be equal to K, that means that K = Pn. Moreover, we can assume that R is small in the sens that P = V p(t) is a good approximation, where V is the volume of R. Hence, under this assumption, one has

p(t) = K t nV , (6.11)
where K t represents the mean number of point falling into R. Then, we can separate R D+2 into a discrete set of bins and count the number of points falling in each bin.

Or we can use a Parzen window implanted at each point and integrate the number of points around each window, which is a smoother estimation. More formally, we have

K hs,x = n i=1 k(t -t i , σ 1 ).
In previous chapter, we used for ABG interpolation a Parzen window method with a Gaussian (or rbf) kernel: k(t it, σ 1 ) = exp( t it 2 /(2.σ 1 )). But, there are several possible kernels that one can consider classicaly in density estimation,

• triangle kernel:

k(u) = (1 -|u|) 1 (|u|≤1) K(u) = (1 -|u|) 1 (|u|≤1) , (6.12)
• Epanechnikov kernel:

k(u) = 3 4 (1 -u 2 ) 1 (|u|≤1) K(u) = 3 4 (1 -u 2 ) 1 (|u|≤1) , (6.13)
• quadratic kernel: • cubic kernel:

k(u) = 15 16 (1 -u 2 ) 2 1 (|u|≤1) K(u) = 15 16 (1 -u 2 ) 2 1 (|u|≤1) , (6.14
k(u) = 35 32 (1 -u 2 ) 3 1 (|u|≤1) K(u) = 35 32 (1 -u 2 ) 3 1 (|u|≤1) , (6.15)
• Gaussian kernel:

k(u) = 1 √ 2π e -1 2 u 2 K(u) = 1 √ 2π e -1 2 u 2 , (6.16)
• circular kernel:

k(u) = π 4 cos π 2 u 1 (|u|≤1) K(u) = π 4 cos π 2 u 1 (|u|≤1) . (6.17)
Let us continue using the Gaussian kernal, such that we can write

K hs,x = n i=1 k(hs -hs i , σ 1 )k(x -x i , σ 1 ).
The multiplication of the two Parzen windows helps us to bring both the spatial and spectral information. Hence the density function is now:

p(t) = p(hs, x) = 1 n n i=1 k(x -x i , σ 1 )k(hs -hs(x i ), σ 1 ), (6.18)
As we seek for P(HS = hs(x)|x, ω), we use Bayes theorem to obtain:

P(HS = hs|x, ω) = n -1 n i=1 k(x -x i , σ 1 )k(hs -hs(x i ), σ 1 ) p(hs, x)dhs . ( 6

.19)

Then by applying the expression (6.8), it is obtained:

f (x) = hs n i=1 k(x -x i , σ 1 )k(hs -hs(x i ), σ 1 ) p(hs, x)dhs dhs = n i=1 n -1 hs k(x -x i , σ 1 )k(hs -hs(x i ), σ 1 )dhs n i=1 k(x -x i , σ 1 )k(hs -hs(x i ), σ 1 )dhs . ( 6

.20)

Using the fact that hs k(hshs(x i ), σ 1 )dhs = hs(x i ), and changing the variable we finally have: .21) In this model, and going back to the problem of multimodal SEM, we do not consider the information brought by the other modalities that can be useful to improve the estimation of hs(x). In order to deal with this idea, and as suggested by [START_REF] Shalizi | Advanced data analysis from an elementary point of view[END_REF] in section 4.3 chapter 4, we modify the single kernel by a product of kernels:

hs(x) = n i=1 hs(x i )k(x -x i , σ 1 ) n i=1 k(x -x i , σ 1 ) . ( 6
k(x -x i ) = j k j (x -x i , σ j ). 144

SPATIAL KERNEL REGRESSION

In particular, the Parzen window kernel will be:

k(x -x i , σ 1 ) = k(x -x i , σ 1 )k(r(x) -r(x i ), σ 1 )k( hs(x) -hs(x i ), σ 1 ), (6.22)
such that this kernel allows us to evaluate the coefficients of the window using the spatial and spectral information, as long as the information of the BSE image r.

The image estimate will be given by:

hs(x j ) = n i=1 hs(x i )ω((x j -x i )σ -1 1 ) (6.23)
with the normalized weights

ω((x j -x i )σ -1 1 ) = k(x j -x i , σ 1 ) n i=1 k(x j -x i , σ 1 )
.

This notation is correct with the rbf kernel.

The question that we can ask ourselves now is what are the statistical properties of such estimator? Indeed, we would like to be sure that if we use another "training

set" N 2 = {x i |i ∈ [1, n]}, with {x i |i ∈ [1, n]} i.i.d. ∼ P (x)
, the regression would stable. Moreover, we do not have an unlimited access to data point and therefore the influence of the number of data should be understood. Finally, we need to know the influence of the (smoothing) scale parameter σ 1 .

Let us denote the real function hs(x) and let hs(x, N 2 ) be the function that is estimated using a Parzen window of parameter σ 1 from training set N 2 . We have the following statistical properties for the normalized kernel function:

w(u)du = 1, uw(u)du = 0, u 2 w(u)du = σ 2 k .
Let us first evaluate the bias of the estimator:

bias 2 = E N 2 hs(x, N 2 ) -hs(x) 2 P (x)dx, (6.24)
thus we have bias

2 = E N 2 n i=1 hs(x i )ω((x i -x)σ -1 1 ) -hs(x) 2 P (x)dx. ( 6

.25)

We can consider that the function hs is locally smooth of first order. This hypothesis is justified by the fact that the function is smoothed by a local filter for denoising it. Under this assumption, a Taylor expansion can be compputed, i.e.,

hs(x

i ) = hs(x) + (x -x i )hs (x) + (x -x i ) 2 2 hs (x) + O((x -x i ) 2 ) (6.26)
Using this approximation in (6.25),

bias 2 = E N 2 n i=1 (x -x i )hs (x)ω((x i -x)/σ 1 ) + (x -x i ) 2 2 hs (x)ω((x i -x)/
E N 2 ( ĥs(x, N 2 ) -hs(x)) = n i=1 (x -x i )hs (x)ω((x i -x)/σ 1 ) + (x -x i ) 2 2 hs (x)ω((x i -x)/σ 1 ) P (x)dx n (x -y)hs (x)ω((y -x)/σ 1 ) + (x -y) 2 2 hs (x)ω((y -x)/σ 1 ) P (y)dy.
By changing the variable z = y-x σ 1 , we have:

E N 2 ( ĥs(x, N 2 ) -hs(x)) = n hs (x) zσ 1 ω(z)P (x + zσ 1 )σ 1 dz + hs (x) zσ 1 ω(z)P (x + zσ 1 )σ 1 dz (6.27)
We also consider that the probability distribution P (x) is locally smooth of second order such that its Taylor approximation:

P (x + zσ 1 ) = P (x) + zσ 1 P (x) + (zσ 1 ) 2 2 P (x) + o((zσ 1 ) 2 ),
makes sense. By injecting the expansion into (6.27) we obtain:

E N 2 ( ĥs(x, N 2 ) -hs(x)) = nσ 3 1 σ 2 k (hs (x)P (x) + 1/2.hs (x)P (x)) + o((zσ 1 ) 2 ). (6.28)
Therefore, we have the final expression for the bias:

bias 2 = (nσ 3 1 σ 2 k ) 2 (hs (x)P (x) + 1/2.hs (x)P (x)) 2 P (x)dx + o((zσ 1 ) 2 ). (6.29) Doing the same kind of calculation for the variance, one has:

variance = E N 2 ĥs(x, N 2 ) -E N 2 ĥs(x, N 2 ) 2 P (x)dx, = σ(x)R(K) nσ 1 P (x) P (x)dx, (6.30)
with R(K) = w(u) 2 du, P (x) is the estimate of P (x) obtained thanks to the Parzen windows, σ(x) represents the noise present on each hs(x) that corrupts the data. See formula (4.11) in [START_REF] Shalizi | Advanced data analysis from an elementary point of view[END_REF] for more details.

Concerning the smoothing parameter, we can see that when its value is small, i.e., σ 1 → 0, so we focus on small scale details and thus we reduce the bias of estimation. In the contrary, the variance of the estimator is increased. Since by reducing σ 1 , with less effect of smoothing, each of our predictions will be an average of less observations, thus they are going to vary more. In summary, this technique has a bias and we just need to find correctly appropiate σ 1 according to the compromise we want to have.

We provide in Figure 6.1 an example of Gaussian kernel regression for a 1D signal. For the case of 2D images, the corresponding local neighbourhood is represented in red in Figure 6.2, thanks to it we approximate the point hs(x) as a weighted mean. More precisely, we can now revisit the expression of the SEM image fusion using the AGB filter, proposed in previous chapter, from the viewpoint of kernel regression:

We can now revisit the expression of the SEM image fusion using the AGB filter, proposed in previous chapter, from the viewpoint of kernel regression:

hs(x) = 1 W p (x) x i ∈E hs(x i ) k(x, x i ), (6.31)
where the AGB kernel is written as the product of two rbf kernels:

k(x, x i ) = k space (x, x i )k EDX+BSE (x, x i ) (6.32)
given by

k space (x, x i ) = g s ( x i -x )M ( O(x i )), k EDX+BSE (x, x i ) = g s (| O(x i ) -O(x)|)g s (|r(x i ) -r(x)|),
with a typical normalization term:

W p (x) = x i ∈E k(x, x i ), (6.33)
M is a mask that allows us to consider just the data on the values of abundance images available at the low resolution:

M ( O(x i )) = 1, if O(x i ) = 0 0, otherwise (6

Gaussian kriging

Gaussian process regression [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF], also called Gaussian kriging, consists in starting with a linear model regression where the value of hs k (x 0 ) is obtained by weighted linear combination of basis functions φ i , 1 i N , i.e.,:

hs(x) = N i=1 ω i φ i (x) + , (6.35)
where the family φ i of N basis functions is used to project the data onto another feature space, and ∼ N (0, σ). By writing φ(x) = (φ 1 (x), . . . , φ N (x)) t ∈ R N and ω = (ω 1 , . . . , ω N ) t , we can put the linear model in vector form as follows:

ĥs(x) = ω t φ(x) + . (6.36)

One can directly see the analogy with the problem of the precedent section, so we have:

P(HS = hs(x)|x, ω) = N (ω t φ(x), σ 2 I) (6.37)
Let us consider a prior over the parameter ω expressing our confidence:

w ∼ N (0, S 0 ), with S 0 ∈ M 1,N (R). We are interested in the value of the joint distribution function of values hs(x 1 ), ..., hs(x n ). From (6.36), we have:

hs(N 2 ) ∼ N (0, φ t S 0 φ + σ 2 I), (6.38)
where hs(N 2 ) = (hs(x 1 ), . . . , hs(x n )) t and φ is a matrix of M n,N (R), such that φ(i, j) = φ j (x i ). We can then define a kernel matrix of the form: K = φ t S 0 φ.

GAUSSIAN KRIGING

Instead of defining basis functions, we can just work with kernels and project the data onto an infinite nonlinear feature space and then apply the representer theorem, see the theory in Chapter 1. We can see from equation (6.38): cov(hs(x i ), hs(x j )) = K(x i , x j ) + σδ(x i , x j ). (6.39) From this expression we can see that nearby data are more correlated, or in other words, the correlation between the data is governed by their distances. This lead to a concept called Gaussian process.

A Gaussian process is a statistical random function f , for which any finite subset of f taken at different spatial locations follows a multivariate Gaussian distribution and it is written as f ∼ GP(0, K).

We consider that hs follows a Gaussian Process, which means that the covariance at two positions satisfies (6.39). Let us assume that HS(x * ) ∼ N (0, K(x * , x * )), which means that there is no noise on it. Moreover the joint distribution of the training outputs hs and the test outputs hs(x * ) is given by

hs(N 2 ) hs(x * ) ∼ N 0, K + σI K * K * K * * (6.40) with K * = φ t S 0 φ(x * ) and K * * = K(x * , x * ).
Then by doing calculation, the posterior probability is obtained as:

P(hs(x * )|x * , hs((N 2 )), (N 2 )) = N (K * (K + σI) -1 hs(x), K * * -K * (K + σI) -1 K t * ). ( 6

.41)

The unbiased estimator for the regression is therefore:

hs(x * ) = E(hs|x * , hs(N 2 ), N 2 ) = K * (K + σI) -1 hs(x). (6.42)
We can see that the prediction is a linear combination of the observations, with a bias: bias = E N 2 hs(x, N 2 )hs(x) 2 P (x)dx. (6.43) We need to integrate:

E N 2 ( hs(x, N 2 )) = E N 2 E(hs|x * , hs(N 2 ), N 2 ) = E(hs|X = x * ) = hs(x * ).
Hence, the solution has no bias. The variance is more complicated to evaluate since we do not have access to the probability of the sampling. That is why, we did not calculate it. We note that Gaussian process regression can be seen as a simple kriging applied in a Gaussian random field. This assumption is quite advantageous when the space coordinate lie on a high dimensional space since it provides an interesting way to estimate the covariance (the kernel). It is a major issue since the Gaussian process regression results depend on the choice of covariance model. In practice, instead of fixing the covariance function, it is commonly assumed to use a parametric family of functions k(x i , x j , θ) = k θ (x i , x j ), and then try to infer θ. A CHAPTER 6. SPATIAL REGRESSION FOR IMAGE PANSHARPENING: APPLICATION TO MULTIMODAL SEM IMAGE FUSION 149 technique used in the community of machine learning to learn the hyperparameters θ are based on maximizing the log likelihood function of p(hs(N 2 )|θ) :

log p(hs(N 2 )|θ, x) = - 1 2 hs(N 2 ) T K -1 θ hs(N 2 ) - 1 2 log det(K θ ) - n 2 log 2π
and maximizing this marginal likelihood towards θ provides the complete specification of the Gaussian process hs. One can briefly note that the first term corresponds to a penalty term for a model's failure to fit observed values and the second term to a penalty term that increases proportionally to a model's complexity. Figure 6.4 provides an illustration of Gaussian process regression for 1D signal. The kernel hyperparameters were learned to optimize the log-marginal likelihood.

Ordinary kriging and pansharpening fusion of information

In the Gaussian process, the model is a linear combination of the observations. The ordinary kriging [START_REF] Matheron | Principles of geostatistics[END_REF] seeks also for this kind of mathematical assumption. It is a model based on a simple linear model regression where the value of hs(x * ) is calculated as a linear combination of the other known realisations:

hs(x * ) = n i=1 ω i (x * )hs(x i ).
However, we do not have access to different realisations of the random field and consequently, this variance. To achieve statistical inference from a single event, the theory of geostatistics has replaced the hypothesis of independent repetitions with some assumption on the random field. The first hypothesis consists in considering that some of its characteristics are identical on all the positions, like the mean for example. Another hypothesis is to consider that the expectations of some quantities are accessible by integrals over space. These assumptions are called stationarity and ergodicity. More precisely, our random field is stationary of second order if we have: E HS (HS(x i )) = E HS (HS), (6.48) cov HS (HS(x i ), HS(x j )) = cov HS (HS(x j + τ ), HS(x j )) = C(τ ). (6.49)

The ergodicity is a property which reinforces the notion of stationary and provides the almost sure convergence of spatial empirical mean to the true mean when the field "goes to infinit".The random field HS is ergodic if one has:

1 |E| E (HS(x)dx) -→ E→+R 2 E HS (HS) (6.50)
where |E| represents the volume of space E, i.e., the space domain of the image.

Thanks to these hypotheses we are able to estimate the covariance of the field. In order to guarantee that this estimator is unbiased, it is necessary to add the following additional condition:

N i=1 w i (x * ) = 1.
Then, we have:

E HS n i=1 w i (x * )HS(x i ) = n i=1
w i (x * )E HS (HS(x i )) = E HS (HS(x i )). (6.51) This leads to a new cost function: We derive it according to w i and µ and equal to zero, and we obtain the following system:    2 n j=1 w j (x * )C(x j , x i ) -2C(x * , x i ) -2µ = 0 ∀i ∈ [i, n] where N τ is the number of elements distant τ . However, working with the empirical variogram does not give satisfactory esults. First, because this function is not necessary negative definite, we might not have a global minimum of the kriging functional. The second point is that we want to have results with a generalization power, thus avoiding an overfitting to the available data. To overcome this issue the basic idea is to fit a model of variogram. There are different admissible parametric variograms,

L(w) = var hs j (x * ) - n i=1 w i (x * )hs(x i ) + 2µ N i=1 w i (x * ) -1 . ( 6
• exponential:

γ(τ ) = C 1 -e -τ a , (6.57)
• Gaussian: In our case, in order to be consistent with the rest of the chappter, we adopted the Gaussian variogram, such that its parameters to be fit are a and C. In addition, there is another parameter which is defined as

γ(τ ) = C 1 -e -(
C 0 = lim τ →0
γ(τ ), called the nugget parameter, and it represents the amount of variance not explained by the model that we have chosen. In order to estimate the parameters, that let us as write as hyperparameter θ, there are different techniques. The simple one that we consider here consists in taking the hyperparameter θ * that minimizes the following functional:

θ * = arg min K k=1 (γ θ (τ k ) -γ(τ k )) 2 , (6.62)
where the experimental variogram is evaluated at K distances τ k . Contrary to the previous section, here we take advantage of the geometrical information, which is crucial in a low dimensional field. Using the experimental variogram as γ θ * (τ ), we can now write the system to perform the kriging estimation: Let us come back now to our problem of multimodal pansharpening and how to inject the information of the BSE image r into the multispectral hs.

       0 γ θ * (
First, we need to impose the assumption that the information is locally stationary of second order. The idea consist in computing a partition of the image domain E into "homogeneous" classes, called in image processing, superpixels. In particular, we have used the SLIC algorithm for superpixels, which is accepted in the stateof-the-art as good algorithm. Then, on each of this superpixel, we perform a local kriging of the abundance maps HS. Let us write SP (x * ) the super pixel containing x * . Hence, our formulation of the local kriging is : Second, the approach to introduce the morphological information of BSE image r into the pansharpening of hs consists in adding a regularization term that would inject this information. A regularization that is often considered on regressionlike is a Bayesian regularization with a Gaussian prior which is linked to the ridge regularization. Here we use a Laplace prior on ω which is related to the Lasso regularization which depends on r. Namely, the distribution of the prior is P(ω|k, λ) = K n j=1 λC r (x * , x j ) exp (-λC r (x * , x j )w j (x * )) , where C r is a covariance function of the BSE image, and K a constant of normalization. Thus, we can rewrite the posterior as: P(ω|hs 1 (x * ), . . . , hs p (x * ), C r , λ, hs(N 2 ), σ) = P(hs 1 (x * ), . . . , hs p (x * )|ω, hs(N 2 ), σ) P(ω|C r , λ) P(hs 1 (x * ), . . . , hs p (x * )) . (6.65)

L(w) = var
The logarithm of the posterior distribution is given by: log(P(w|Z 0 , Z, σ, k, λ) C r (x * , x j )w j (x * ), (6.67)

We need again the convex constrain: where µ is a Lagrange parameter. For the BSE image, we also use Gaussian model for the covariance:

C r (h) = exp - h a ,
where a should be fit from the empirical covariance. The parameter λ is fundamental in this (multimodal) regularized kriging since it represents the quantity of spatial information of r to be injected into hs.

6.5 Results and discussion

Evaluation criteria of pansharpening algorithms

The criteria used to evaluate the quality of the merged information are those conventionally considered in the context of the pansharpening literature. We will use 156 6.5. RESULTS AND DISCUSSION exactly the same criterion that we developed in the previous chapter. More precisely, we considered the following criteria: C1: The spectral distortion between the enhanced multispectral image and the real multispectral image at the nominal resolution should be as small as possible. Or, in other terms, we would like to find the same materials for each pixel as the original image at high resolution;

C2: The spatial distortion between the enhanced multispectral image and the real one should not be too high.

C3: The injected information. Since the information provided by the various modalities can be relatively different, we would like to inject the useful one to improve segmentation and characterisation of the EDX image.

Let us write HS the multispectral abundance EDS image at the same resolution that R, that has been provided by the sensor at high resolution. In a way it is the ground truth such that our enhanced images HS should be compared to HS. Moreover, we denote by M ∈ M n,D (R) and M ∈ M n,D (R) the two matrices representing respectively HS and HS, where n is the total number of pixels (i.e., n = N 1 × N 2 ), and D the number of abundance maps. We write by M i,: , ∀i ∈ [1, n], a spectra of HS, and M :,k ∀k ∈ [1, D] a map of HS. We have now all the notations to introduce the five parameters.

Cross correlation (CC) is a measure that evaluates the spatial distortion defined as:

CC( HS, HS) = 1 D D k=1
CCS( M :,k , M :,k ), (6.71) where:

CCS( M :,k , M :,k ) = ( n i=1 M i,k -µ M :,k )( n i=1 M i,k -µ M :,k ) n i=1 ( M i,k -µ M :,k ) 2 n i=1 ( M i,k -µ M :,k ) 2
, with µ M :,k = n -1 n i=1 M i,k being the empirical mean. The CC is optimal when it is close to 1. where < •, • > is the Euclidean inner product of vectors associated to the norm L 2 , and where • is the L 2 norm of vectors. The SAM is optimal when it is near to 0.

Spectral Angle

Root mean squared error (RMSE) measures the mean residual error of fusion and is defined as: RM SE( HS, HS) = M -M F n.D , (6.73) where • F is the Frobenius norm of a matrix A, i.e., A F = trace(AA t ). The RMSE is optimal when it is near to zero. The ERGAS is optimal when close to 0.

Cross correlation gradient (CC g ) is a measure that evaluates the spatial distortion defined as where Rg represents the image morphological gradient [START_REF] Soille | Morphological image analysis: Principles and applications[END_REF] of the BSE image converted into a vector and Mg :,k is the image morphological gradient of the enhanced abundance map at the nominal scale converted also into a vector.

An optimal enhancement method would be a compromise between the criterion C3 and the C1 and C2.

Evaluation on dataset 1

The EDS abundance maps of dataset 1 are provided in Figure 6.5. The results of the enhanced abundances HS obtained by the different techniques provided in Figures ??, 6.7 6.13 and 6.9, where the abundances are visualized as RGB color images. Quantitative results for this case according to the five measures are presented in Table 6.1.

From these results, we note that AGB presents good performance but it smoothes a lots the image, while kriging techniques inject more gradually geometric/textural information from the BSE image. Then, even if based on these results it is difficult to definitely compare these techniques, we might see that the advantage of the kriging technique is it flexibility since the operator just need to select the quantity of parsimony based on the backscattered he wants to apply.

Evaluation on dataset 2

In the case of this multimodal SEM dataset, which has a higher level of noise, we try to inject the BSE image R, depicted in Figure 6.10, to increase the resolution of the EDS abundance maps provided in Figure 6.11, and also to denoise the images. The results for our ABG method are given in Figure 5.14. For comparison, Figures 6.12 and 6.13 provide a visualization of the results obtained using the other methods. Quantitative results are given in Table 6.2, which lead to similar conclusions as for dataset 1.

Techniques

Conclusion

In this chapter we reformulate the pansharpening problem in the form of a spatial regression problem. We have proposed several solutions to deal with such regression. We compared and evaluated these different solutions for the specific problem of enhanced EDS in multimodal SEM imaging. An interesting point of these two techniques is that they work on the full vector space, such that the full spectra of the multispectral image are used. The kernel regression technique preserves the positivity of the data too. The latter is physically important since the data we are dealing with represents quantities of materials. The kriging technique provides results that are visually interesting. It appears that the technique based on regularized kriging seems to be easier to use since there is just a single parameter to tune that corresponds to the quantity of the high spatial resolution information to be injected. The technique based on kernel regression seems to have better results on the test datasets, with the criteria that we used. However this technique have four parameters, which are on the one hand, the size of the kernels, and on the other hand the three scaling parameters σ of each rbf kernel. In order to improve this technique, it might be interesting to try to learn some of these parameters automatically. Chapter 7

Conclusions and Perspectives

Summary of main contributions

In this thesis we addressed the problem of finding relevant representations for multivariate and multimodal images.

The idea that we developed for this kind of representations is to include spatial information in the statistical learning process. We have in particular worked on multivariate spatial dimensionality reduction, spatial classification and multimodal spatial regression.

Hence, the main contributions of this thesis are just based on an appropiate use of spatial information in machine learning techniques in order to fight the curse of dimensionality on the manifold of multivariate image pixels.

We have use methods and algorithms from mathematical morphology, wavelet and scattering transform, geostatistics, regression and kernel-based machine learning. We compared the proposed techniques theoretically and also in terms of experimental results with state-of-the-art methods.

To summarize, the contributions of our work include the following points.

Chapter 3. We have shown how to deal with spatial information in the process of multi/hyper-spectral image dimensionality reduction thanks to mathematical morphology operators. The representation techniques that we introduced in this chapter are based on the notion of the MorphPCA. Dimensionality reduction is done on the space of the data without any approximation regarding the notion of morphological covariance V Morpho , where V Morpho is a covariance handling the morphological/spatial relationships between the image bands of the multi/hyper-spectral image. We proposed different alternatives to compute these covariance matrices. Globaly, the technique is simple in terms of computation and memory storage. Moreover, we also proposed some criteria to assess the quality of the image representation after the dimensionality reduction. Some of them are based on mathematical morphology, namely the 3D pattern spectrum and the α-flat zone partition, and are used to evaluate if the reconstructed image preserves globally and locally the similarity to the original hyperspectral image. Finally, we also consider the interest of Mor-phPCA for supervised classification on the reduced data. According to the entire set of criteria, adding spatial information improves the dimensionality reduction. However, as we have shown, a good dimensionality reduction is 7.2. SUGGESTIONS FOR FUTURE WORK obtained when combining spatial and spectral feature spaces.

Chapter 4. We have addressed the interest of spatial information in the process of multi/hyper-spectral image classification. To perform this classification of pixels, we implemented a deep convolutional descriptor. This descriptor is based on two innovative techniques called scattering transform and kernel mean map. Thanks to it, a translation-invariant representation of the texture is obtained, which is also Lipschitz stable to deformations. To evaluate this descriptor, we just calculated the mean of the scattering coefficients provided by each layer embedded on the Hilbert space of a rbf kernel. This classifier simplifies convolutional neural network since fewer parameters are needed. This descriptor constitutes an interesting way to perform deep learning when the training data set is limited. We studied the statistical and pattern recognition properties of this descriptor, and prove also empirically that it presents a good classification performance.

Chapters 6 and 7. We have studied how spatial information is naturally used on a process of fusion of information called pansharpening. To perform this fusion of information of multimodal SEM images, different techniques were deeply considered. One of the major difficulties of this part was the lack of literature on this kind of techniques applied to multimodal SEM. Thus, our starting point was the pansharpening state-of-the-art solutions proposed for multimodal remote sensing images. We have developed three techniques. The first one consists in using the nonparametric kernel regression, where the kernel handles spatial and spectral information. One of these well known kernels widely used in image processing is the bilateral filter. This approach presents a good compromise according to the criteria of evaluation that we used. We proposed also to use a technique based on morphological wavelets for upsampling. Finally, the third approach is based on ordinary kriging. From our viewpoint, the last is quite interesting since we have a way to measure how much of the second modality we want to inject into the first one. However, the inconvenience with such a technique is that we lose some physical properties of the final results, even if visually the results seem better.

Suggestions for future work

We finally discuss some open questions and suggestions based on the works we proposed.

Chapter 3: The following points deserve some developments,

• The MorphPCA dimensionality reduction techniques proposed were used on hyperpectral images. It might be interesting to consider their interest for other kind of multivariate images, including also multimodal ones.

• Dimensionality reduction is used to denoise images sequences, and thus investigating the interest of MorphPCA approaches might be interesting in time series.

• Some fusion of information techniques are based on the PCA, MorphPCA could replace the PCA on these techniques.

• It might be interesting to theoretically link MorphPCA with kriging, in order to have a geostatistical theory of spatial dimensionality reduction.

Chapter 4: The potential of this part of the thesis is important and different aspects should be investigated in ongoing work.

• To use other classification techniques than the SVM. We would suggest to apply Gaussian process, since it is linked with kriging and with the scattering.

• To focus more on learning the perfect Hilbert space where the scattering transform is embedded.

• To study the interest of this descriptor for the simulation of hyper/multi-spectral textures.

• To use this descriptor to solve problems other than classification.

• In the process of estimating the rbf feature space, we need to project the data onto a random space. Estimating this space in a non-random way might improve the results. By doing that, we might not emmbed the data on the rbf feature space and we can expect an improvement of the precision.

• Another interesting solution might be to use these descriptors as local visual descriptors (replacing SIFT for example), and use them in typical computer vision tasks.

• Instead of using a Gaussian average, one could use a perceptron algorithm, and learn the coefficient of the weights, closer to the convolutional neural network paradigm.

Chapters 6 and 7. Related to the problem of image fusion and pansharpening, we suggest for future investigations to work on:

• To assess the interest of using other SEM modalities, e.g., the secondary electron image, in the process of image fusion;

• To study how to improve the selection of the EDX pixels images carefully, so that when we want to increase the dimension, the most important information is used;

• To study if there exists a link between spatial regression and with compressive sensing theory;

• To develop a kriging model with an additional constraint of nonnegative coefficients (nonnegative kriging).

Résumé

Mots Clés

Figure 1 :

 1 Figure 1: A set of 900 data points of dimension D = 2. In (a) the data, in (b) their d i,j distance matrix, in (c) the k = 3 clusters of the dataset obtained by k-means.

Figure 2 :

 2 Figure 2: A set of 900 data points of dimension D = 100. In (a) the d i,j distance matrix of the dataset. As we can see, it is more difficult to separate some classes. In (b), the k = 3 clusters of the data obtained by k-means: clustering is not correct because of the curse of dimensionality.

Figure 3 :

 3 Figure 3: In this figure we have selected randomly 500 data points,v i ∈ [0, 1] D , i = 1, • • • , 500, where D is the dimension of the feature space. We plot R =

Figure 1 . 1 :

 11 Figure 1.1: Hyperspectral image [98]

CHAPTER 1 .Figure 1 . 5 :

 115 Figure 1.5: (a) Synthetic manifold with two concentric spheres, (b) polynomial KPCA with p = 5, (c) Gaussian KPCA with σ, (d) Gaussian KPCA with 5σ,(e) Gaussian KPCA with 8σ, (f) Gaussian KPCA with 15σ.

Figure 1 . 6 :

 16 Figure 1.6: (a) Flower-like synthetic manifold, (b) Polynomial KPCA with p = 5, (c) Gaussian KPCA with σ, (d) Gaussian KPCA with 5σ,(e) Gaussian KPCA with 8σ, (f) Gaussian KPCA with 100σ.

35 Figure 1 . 7 :

 3517 Figure 1.7: Illustration of the geometry of the hyperplane for binary classification.

Figure 1 . 9 :

 19 Figure 1.9: Illustration of the neural network with one input layer, one output layer and one hidden layer.

Figure 1 . 11 :

 111 Figure 1.11: Illustration of a simple Convolutional Neural Network layer [16].

Figure 1 . 12 :

 112 Figure 1.12: Illustration of a typical Convolutional Neural Network architecture [2].

Figure 2 Figure 2

 22 Figure 2.2: (a) Channel number 50 of Pavia hyperspectral image and (b) its morphological decomposition by area openings γ a s l , s l = {0.5 10 2 , 1 10 2 , 5 10 2 , 7 10 2 , 1 10 3 , 2 10 3 , 5 10 3 , 7 10 3 , 1 10 4 , 1.2 10 4 , 1.5 10 4 , 2.5 10 4 }. Last image in (b) corresponds to γ a s S , s S = 2.5 10 4 ; the other images in (b) are (γ a s l-1 (f )γ a s l (f )). Note that the contrast of images has been enhanced to improve visualization.

Figure 2. 4 .

 4 We propose three variants of MorphPCA which are summarized in the flowchart of Figures 2.5, 2.6 and 2.8. An example of three different bands embedded in the space produced by these MorphPCA techniques is depicted in Figure2.9.

Figure 2 . 4 :

 24 Figure 2.4: Global process of MorphPCA.

Figure 2 . 5 :

 25 Figure 2.5: Process of scale-space decomposition MorphPCA.

53 Figure 2 . 6 :

 5326 Figure 2.6: Process of pattern spectrum MorphPCA.

Figure 2 54 2Figure 2 . 8 :

 25428 Figure 2.7: (a) Example of pair of binary images for pattern spectrum correlation discussion. (b) Example of triplet of binary images for distance function correlation discussion.

Figure 2 . 9 :Figure 2 . 10 :

 29210 Figure 2.9: Top, three examples of spectral bands of Pavia image: (a) 1, (b) 50, (c) 100; middle, (d), (e), (f) PS of corresponding spectral bands; (g), (h), (i) Molchanov distance functions of corresponding spectral bands.

Figure 2

 2 Figure 2.11: (a) A 3-variate image (first three eigenimages after PCA on Pavia hyperspectral image) and (b) its corresponding α-flat zone partition into 84931 spatial classes using the Euclidean distance.

Figure 2 . 12 :

 212 Figure 2.12: RGB false color visualization of first three eigenimages from Pavia hyperspectral image: (a) classical PCA on spectral bands, (b) scale-decomposition MorphPCA, (c) pattern spectrum MorphPCA, (d) distance function MorphPCA.

Figure 2 .

 2 Figure 2.13: Hyperspectral band projection into the first two eigenvectors (i.e., image manifold) from Pavia hyperspectral image: (a) classical PCA on spectral bands, (b) scale-decomposition MorphPCA, (c) pattern spectrum MorphPCA, (d) distance function MorphPCA.

  (a) for Pavia hyperspectral image, (b) for Houston hyperspectral image, (d) for Indian Pines hyperspectral image. The values have been normalized to the worst case, which gives 100.

Figure 2 . 14 :

 214 Figure 2.14: Intrusion/Extrusion parameters for PCA and the different variants of MorphPCA from Pavia hyperspectral image: (a) Q(K), (b) B(K).

Figure 2 . 15 :

 215 Figure 2.15: (a) 3D pattern spectrum distribution of Pavia hyperspectral image and of the different reduced images into d components. (b) Corresponding 3D cumulative pattern spectrum distributions.

  Ground truth PCA, OA : 43.9 Morpho-1, OA : 50.5 Morpho-2, OA : 41.5 Morpho-3, OA : 51.3 Morpho-4 β=0.2 , OA : 43.5 Morpho-4 β=0.5 , OA : 43.1 Morpho-4 β=0.8 , OA : 43

Figure 2 . 16 :

 216 Figure 2.16: Results of supervised classification using least square SVM with a linear kernel on Indian Pines hyperspectral image. Note the OA is the overall accuracy.

Figure 2 . 17 :

 217 Figure 2.17: Results of supervised classification using least square SVM with a rbf kernel on Indian Pines hyperspectral image. Note the OA is the overall accuracy.

Figure 2 . 19 :

 219 Figure 2.19: Results of kappa statistic for the least square SVM with a rbf kernel and different percentage of training set on Indian Pines hyperspectral image, the dimension of the reduced space is equal to 5.

  Figure 3.1: Hyperspectral patch P (s) v (f ).

Figure 3 Figure 3

 33 Figure 3.2: (a) The ground truth for the "Indian Pine" dataset. (b) The local standard variation for this image, which represents the local spatial variability of the spectral image.

Figure 3. 4 :

 4 Figure 3.4: The different steps of the deep weigthed mean map scattering process .

Figure 3 . 5 :

 35 Figure 3.5: Process of scattering transform on a signal.

  .6.

Figure 3 . 6 :

 36 Figure 3.6: The deep weighted mean map process.

CHAPTER 3 .Figure 3 . 7 :

 337 Figure 3.7: (a) The difference of two characteristic functions Φ Px and Φ Py . (b) A possible choice of Λ(ω) for a possible kernel. (c) As we can see, this kernel captures the difference between the two characteristic functions.

  Figure 3.8: A representation of the SVM kernel embedding of the distributions, based on a finite data set partitioned into bag of distributions.

  our training set and classification algorithm. So to improve it, we split the training set into two sets. One is called the training set, and the other one the validation set. Then we learn a model on this new training set and then classify the validation set testing a family of parametersσ {0.01 × σ, 0.1 × σ, 0.5 × σ, σ, 3 × σ, 6 × σ, . . . , 30 × σ}.After having performed a sufficient number of validation classifications, we choose the parameter that brings the best results. The issue with this technique is that the result of the choice of the parameter depends on the initial training set. In the case of remote sensing, most of the time the training set is small. To overcome this issue, we propose to build a new kernel mean map by combining a predefined set of kernels mean maps. Let us consider that we have a set of M kernels K = {k 1 , . . . , k M }. One can consider as a multiple kernel the sum kernel : K = M m=1 k m ; or the weighted sum kernel : K = M m=1 β m k m where the coefficients β m can be fixed or learned during the classification process. This leads to multiple kernel supervised learning techniques.

Figure 3 . 9 :

 39 Figure 3.9: In (A) the zeros and first order scattering before the mean, U of an abondance map are represented. On the first line the abondance map is given, and then each line is composed of the 8 orientations at a given scale. In (B) we represent the second order scattering before the mean, U . In (C) and (D) we represent respectively the weigthed deep mean map first layer and second layer.

Figure 3 . 10 :

 310 Figure 3.10: Classification maps for the Indian Pines hyperspectral image using different approaches, with just 5 points per class in the training set. In (a) ground truth, (b) the linear SVM, (c)kernel SVM with Morphological Profile, (d) kernel SVM with K M P × Kmm and s = 15, (e) Scattering transform , (f) deep mean map Scattering transform S H m = 1 and s = 10.

Figure 3 . 11 :

 311 Figure 3.11: Classification maps for the Pavia hyperspectral image using different approaches, with just the classical training set. In (a) ground truth, (b) the linear SVM, (c)the estimated RBF SVM, (d) kernel SVM with KM M and s = 15, (e) kernel SVM with KCN and s = 13.

103

 103

104 4 .

 4 Figure 4.1: Scanning Electron Microscope overview [3].

CHAPTER 4 . 105 Figure 4

 41054 Figure 4.2: SEM interations [117].

Figure 4 . 3 :

 43 Figure 4.3: Backscattered electrons sensors [1].

Figure 4 . 4 :

 44 Figure 4.4: Emission domains of various interactions[117].

  Figure 5.1: Multimodal SEM image acquisition. In the red square ( a,b,c,d) we represented 4 abundances. In the blue square (e) we represented the backscattered image and in the green square (f) the secondary electron image.

  Backscattered electrons were also registered to obtain a simulated BSE image. The image simulations were performed with 10000 electrons trajectories for each pixel with varying image size and forcing factor F . Varying F allowed to tune the intensity of Poisson noise of the simulated images. Two images of 1024 × 1024 pixels with F = 1 and F = 100 and three images of 256 × 256 pixels with F = 1, F = 10 and F = 100 were simulated. Typical simulation time was about 10 days for a 1024 × 1024 pixels image on an 8 processors working station.Figure 5.2 depicts the ground truth sample used in the simulation: it is composed of regular geometric shapes of pure and binary composition with low mean atomic number element (Al 2 O 3 , SiO 2 ), medium (Fe, Co, Ni) and high (Pt).

Figure 5 .

 5 3 provides the simulated total backscattered electrons image. It is worth noticing that phases with close mean atomic number (Al 2 O 3 and SiO 2 ; Fe, Co and Ni) are hard to be distinguished in the BSE image. In Figure5.4(a) is given the pixelwise mean image of the simulated EDX map. We can notice that the pixelwise mean spectrum reveals clearly the Pt zone due to high Bremsstrahlung emission in this high mean atomic number zone.

Figure 5 .

 5 4(b) shows the image corresponding to the Al Kα energy channel, the one having the less Poisson noise.

Figure 5 . 2 :

 52 Figure 5.2: Ground truth of the simulated multimodal SEM image.

Figure 5 .

 5 Figure 5.3: Simulated backscattered electron image.

CHAPTER 5 .Figure 5 . 4 :

 554 Figure 5.4: (a) The pixelwise mean image of the simulated EDX image (256 × 256, F = 1). (b) Simulated Al Kα energy channel.

Figure 5 . 5 :

 55 Figure 5.5: Model of a hierarchical decomposition of the information.

118 5. 4 .

 1184 FUSION OF SEM INFORMATION BY ABUNDANCE GUIDED BILATERAL FILTER (AGB)

Figure 5 . 6 :

 56 Figure 5.6: The two images used for the fusion of SEM information by AGB filter.

CHAPTER 5 .

 5 ENHANCED EDX IMAGES BY FUSION OF MULTIMODAL SEM IMAGES USING PANSHARPENING TECHNIQUES 121

5 . 25 )

 525 The adjoint erosion of the dilation sampling operator is given by σK (f S ) :fun(F S ) -→ fun(F ) σK (f S )(s) = inf{f S (s)|s ∈ K * (x)} (5.26)and the adjoint dilation of the erosion sampling operator as νK (f S ) : fun(F S ) -→ fun(F ) νK (f S )(s) = sup{f S (s)|s ∈ K * (x)} (5.27)

29 )CHAPTER 5 .

 295 ENHANCED EDX IMAGES BY FUSION OF MULTIMODAL SEM IMAGES USING PANSHARPENING TECHNIQUES 123

  Mapper (SAM) is a measure that assesses the spectral distortion by computing SAM ( HS, HS) = 1 n n k=1 SAM ( M i,: , M i,: ), (5.36) with SAM ( M i,: , M i,: ) = arccos   < M i,: , M i,: > M i,: M i,:   , CHAPTER 5. ENHANCED EDX IMAGES BY FUSION OF MULTIMODAL SEM IMAGES USING PANSHARPENING TECHNIQUES 125

126 5. 6 .

 1266 RESULTS AND DISCUSSION

Figure 5 . 9 :Figure 5 . 10 :Figure 5 . 11 :Figure 5 .Figure 5 . 14 :Figure 5 . 15 :Figure 5 . 16 :

 595105115514515516 Figure 5.9: Zoom of the RGB color image from abundances 1,2,4 of the dataset 1 for the different pansharpening techniques.

CHAPTER 5 .Figure 5 . 17 :

 5517 Figure 5.17: Extraction of the most abundant material at each pixel for the different pansharpening techniques.

Figure 6 . 1 :

 61 Figure 6.1: Illustration of the signal estimation thanks to the Gaussian kernel regression approach. In (a), one can see in black the target function, and in blue the availble data or the training set, which is composed of 30 points. In (b), different training sets of 29 points have been generated based on one of 30 and we perform a regression in the space. Dark curve is the mean prediction over all the different training sets, and the grey are is the variance.

. 34 )

 34 CHAPTER 6. SPATIAL REGRESSION FOR IMAGE PANSHARPENING: APPLICATION TO MULTIMODAL SEM IMAGE FUSION 147 The function O was introduced in the previous chapter and represents the level of mixing of the different abundances. Finally, g s is a Gaussian function of scale parameter s.

Figure 6 .

 6 Figure 6.2: Spatial kernel regression with image information.

Figure 6 . 3 :

 63 Figure 6.3: Illustration of the signal estimation thanks to the Gaussian process regression. In (a), one can see in black the target function, and in blue the availble data or the training set, which is composed of 30 points. In (b), different training sets of 29 points have been generated based on one of 30 and we perform a regression in the space. Dark curve is the mean prediction over all the different training sets, and the grey are is the variance.

  x * )w i (x * )C(x j , x i ) x * )C(x * , x i ) + C(x * , x * ) + 2µ N i=1 w i (x * ) -1 .

N

  j=1 w j (x * ) = 1(6.53) Instead of solving the kriging with the covariance, one can therefore do it with the variogram. The advantage of working with the variogram is we do not lean on the estimation of the mean. Morerover one can define the empirical variogram by: i ) -HS(x i + τ ) 2 2 , (6.56)

Figure 6 . 4 :

 64 Figure 6.4: Illustration of the signal estimation thanks to the kriging regression. In (a), one can see in black the target function, and in blue the availble data or the training set, which is composed of 30 points. In (b), different training sets of 29 points have been generated based on one of 30 and we perform a regression in the space. Dark curve is the mean prediction over all the different training sets, and the grey are is the variance.

N

  j=1 w j (x * ) = 1. Therefore, the use of Lagrange multiplier theorem provides us:x * )w i (x * )C(x j , x i ) -2 n i=1 w i (x * )C(x * , x i ) +C(x * , x * ) -2µ( N i=1 w i (x * ) -1) -2λ n i=1 C r (x * , x i )w i (x * )which, deriving as previously according to w i and µ, involves the following kriging system:   n j=1 w j (x * )C(x j , x i ) -2C(x * , x i ) -2µ -2λC r (x * , x i ) = 0 ∀i ∈ [i, n] 1x 1 ) . . . C(x 1x n 1x n ) . . . C(x nx n 1x * ) -λC r (x * , x 1 ). . . C(x 1x * ) -λC r (x * , x 1 ) 1

  Mapper (SAM) is a measure that assesses the spectral distortion by computing SAM ( HS, HS) = 1 n n k=1 SAM ( M i,: , M i,: ), (6.72) with SAM ( M i,: , M i,: ) = arccos   < M i,: , M i,: > M i,: M i,:   , CHAPTER 6. SPATIAL REGRESSION FOR IMAGE PANSHARPENING: APPLICATION TO MULTIMODAL SEM IMAGE FUSION 157

 2 ,

 2 Synthetic adimensional global error (ERGAS) offers a global measure of the quality of an enhanced image. It is given by the expression:ERGAS( HS, HS) = 100d D k=1 RM SE( M :,k , M :,k ) µ M :,k  where d is the ratio between the linear resolution of the BSE image R and the abundances EDS image HS, i.e., d = R-linear spatial resolution HS-linear spatial resolution .

  , Mg :,k ),

Figure 6 . 5 :Figure 6 . 6 :Figure 6 . 7 :Figure 6 . 8 :

 65666768 Figure 6.5: The six EDS abundance maps of SEM dataset 1 at the nominal resolution.

Figure 6 . 9 : 163 Techniques

 69163 Figure 6.9: Zoom of the RGB color image from abundances 3,5,6 of the dataset 1 for the different spatial interpolation techniques.

Figure 6 . 10 :

 610 Figure 6.10: BSE image of SEM dataset 2.

Figure 6 . 11 :

 611 Figure 6.11: Four EDS abundance maps of SEM dataset 2 at the nominal resolution.

CHAPTER 6 .Figure 6 . 12 :

 6612 Figure 6.12: RGB color image from abundances 2,3,4 of the dataset 2 for the different spatial interpolation techniques.

Figure 6 .

 6 Figure 6.13: Zoom of the RGB color image from abundances 2,3,4 of the dataset 2 for the different spatial interpolation techniques.

  INTRODUCTION de données tridimensionnelles, où chaque image représente une bande spectrale. En général, le nombre de bandes est environ dix pour des images multispectrales, une centaine pour des images hyperspectrales. Il est donc nécessaire d'utiliser des techniques d'apprentissage statistique pour mieux comprendre ces données. Dans ce chapitre, nous nous concentrerons sur un examen des techniques d'apprentissage statistique qui ont été élaborées, afin d'une part, de réduire la dimension de nos données et, d'autre part, de classifier. De nombreuses techniques de réduction de la dimension ont été proposées, nous nous concentrons ici sur l'analyse des composantes principales (APC), multidimensionnelle scalling (MDS), Isomap, Kernel-PCA, diffusion map, Locally Linear Embedding (LLE), Laplacian Eigenmaps and Hessian LLE. Nous rappelons également certains algorithmes d'apprentissage de dictionnaire. Enfin, nous expliquons la régression linéaire, la classification linéare, le Support Vector Machine (SVM) et les Convolutional Neural Networks (CNN) pour la classification.

Ce type de techniques est très utile dans le traitement d'images hyperspectrales. Ces images proviennent de capteurs hyperspectrales et peuvent être considérées comme un cube 11 1.1.

Table 2 .

 2 2: Comparison of PCA and MorphPCA analysis using criteria C1 and C2:

Table 2 .

 2 3: Comparison of hyperspectral supervised classification on PCA and Mor-phPCA spaces using least square SVM algorithm and different kernels: (a) Pavia hyperspectral image, (b) Indian Pines hyperspectral image.

	Overall Accuracy Overall Accuracy Kappa statistic
	with linear kernel with rbf kernel	with rbf kernel
	V V Morpho-1 V Morpho-2 V Morpho-3 V Morpho-4 β , β = 0.2 61.0 ± 1.73 51.51 ± 0.9 59.6 ± 2.2 56.99 ± 1.1 59.9 ± 2.5 V Morpho-4 β , β = 0.5 59.9 ± 1.5 V Morpho-4 β , β = 0.8 57.87 ± 3	84.9 ± 3.1 85.8 ± 2.6 85.2 ± 2.1 86.0 ± 1.9 85.2 ± 1.1 84.6 ± 1.0 84.7 ± 2.5 (a) Pavia image	0.84 ± 1 × 10 -4 0.84 ± 1 × 10 -4 0.84 ± 1 × 10 -4 0.84 ± 1 × 10 -4 0.83 ± 1 × 10 -4 0.83 ± 1 × 10 -4 0.83 ± 2 × 10 -4
	Overall Accuracy Overall Accuracy Kappa statistic
	with linear kernel with rbf kernel	with rbf kernel
	V V Morpho-1 V Morpho-2 V Morpho-3 V Morpho-4 β , β = 0.2 43.5 ± 3.3 43.9 ± 3.6 50.5 ± 3.8 41.5 ± 3.8 51.3 ± 3.2 V Morpho-4 β , β = 0.5 43.1 ± 2.9 V Morpho-4 β , β = 0.8 43.0 ± 2.2 (b) Indian Pine image 75.2 ± 3.7 79.6 ± 3.7 66.6 ± 4.6 79.1 ± 3.2 75.1 ± 2.3 71.2 ± 2.6 69.7 ± 3.3	0.73 ± 4.3 × 10 -4 0.78 ± 4 × 10 -4 0.63 ± 4.5 × 10 -4 0.77 ± 3.7 × 10 -4 0.72 ± 2.6 × 10 -4 0.68 ± 3 × 10 -4 0.67 ± 3.9 × 10 -4

Table 3 .

 3 1: Overall accuracy (OA), kappa statistic, and average accuracy (AA) obtained for different kernels, applied on the AVIRIS Indian Pines hyperspectral data set. 20 Monte Carlo simulations were run. The training set is just 5 samples per class.

Table 3 .

 3 2: Overall accuracy (OA), kappa statistic, and average accuracy (AA) obtained for different kernels, applied on the AVIRIS Indian Pines hyperspectral data set. 20 Monte Carlo simulations were run. The training set is just 15 samples per class.

	kernel	parameters OA	kappa statistic AA
	Linear kernel rbf kernel Morphological profile K M P × Kmm Scattering transform U m = 1 Scattering transform U m = 2 Scattering transform S m = 1 Scattering transform S m = 2 Deep mean map Scattering transform S H m = 1 s = 10 s=15 s = 5 s = 5	75.4 ± 1.1 77.3 ± 1.0 85.1 ± 2.2 92.8 ± 0.5 90.27 ± 0.8 88.60 ± 1.0 72.3 ± 1.2 74.4 ± 1.1 83.1 ± 2.4 91.1 ± 0.6 91.0 ± 1.0 89.1 ± 1.1 91.1 ± 0.9 90.2 ± 1.1 92.6 ± 0.7 91.7 ± 0.9 96.2 ± 0.5 95.9 ± 0.6	76.1 ± 1.9 77.1 ± 1.9 85.6 ± 1.5 93.0 ± 0.6 90.2 ± 0.6 91.7 ± 0.6 92.4 ± 0.7 93.8 ± 0.5 96.08 ± 0.4

Table 3

 3

	kernel	parameters OA	kappa statistic AA
	Linear kernel		73.2	66.6	78.5
	rbf kernel		75.1	67.6	82.6
	Morphological profile		97.1	96.2	96.7
	K M P × Kmm CNN	s=15	97.4 ± 0.6 96.4 ± 0.7 76 74	97.3 ± 0.6 73
	Scattering transform U m = 2 Scattering transform S m = 2 Scattering transform S m = 2 Deep mean map Scattering transform S H m = 1 s = 10 s = 3 s = 15	92.27 ± 3.1 89.60 ± 2.8 95.23 ± 3.3 93.63 ± 3.2 93.53 ± 3.1 91.34 ± 3.4 96.5 ± 2.7 95.7 ± 2.8	94.69 ± 3.3 96.38 ± 3.3 93.18 ± 3.1 96.2 ± 2.6

.3: Overall accuracy (OA), kappa statistic, and average accuracy (AA) obtained for different kernels, applied on the AVIRIS Indian Pines hyperspectral data set. 20 Monte Carlo simulations were run. The training set is just 50 samples per class.

Table 3 .

 3

4: Overall accuracy (OA), kappa statistic, and average accuracy (AA) obtained for different kernels, applied on the University of Pavia hyperspectral data set. The predefined training set was used.

Table 3 .

 3

	kernel	parameters OA	kappa statistic AA
	Linear kernel Scattering transform U m = 1 Deep mean map scattering transform S H m = 1 s=3 γ = 2	62.9 ± 3.8 57.6 ± 4.1 65.1 ± 3.9 60.6 ± 4.2 80.0 ± 4.1 73.7 ± 4.3	51.2 ± 3.9 55.1 ± 3.7 72.2 ± 4.1

5: Overall accuracy (OA), kappa statistic, and average accuracy (AA)obtained for different kernels, applied on the Zurich data set. 20 Monte Carlo simulations were run. The training set is composed of 19 images and the testing set of 1 image. set. Then, on the other 19 images, we perform a validation by selecting randomly 3 images on the validation set and 16 on the training set such that on the training set we have the 8 classes. Once the parameters are learned, we train on the 19 images and then classify on the testing set. The results of the different algorithms are given on Table

3

.5. As one may see, our descriptor produces relevant results in comparison with the others.

Table 5 .

 5 1: Comparison of pansharpening algorithms for enhanced EDS abundance images of dataset 1.

  Figure 5.8: RGB color image from abundances 1,2,4 of the dataset 1 for the different pansharpening techniques. CHAPTER 5. ENHANCED EDX IMAGES BY FUSION OF MULTIMODAL SEM IMAGES USING PANSHARPENING TECHNIQUES 129

	CHAPTER 5. ENHANCED EDX IMAGES BY FUSION OF MULTIMODAL 128 5.6. RESULTS AND DISCUSSION
	SEM IMAGES USING PANSHARPENING TECHNIQUES		127
		Reference Reference		Input image	Input image
		abundance 1			abundance 2
		GS method GS method		PCA method PCA method
		abundance 3			abundance 4
		SFIM method SFIM method		MTF GLP method MTF GLP method
		Techniques CC SAM RMSE ERGAS CC g
		GS	0.23 16.3	1.6	20.3	0.66
		PCA SFIM abundance 5	0.17 17.4 0.14 15.9	1.6 1.74	20.6 27.9 abundance 6	0.67 0.17
		MTF GLP 0.30 15.8	1.4	20.1	0.61
	GFPCA Figure 5.7: The six EDS abundance maps of SEM dataset 1 at the nominal resolu-0.26 15.9 1.5 20.4 0.31
	tion.	AGB	0.30 15.7	1.4	20.1	0.45
		BGM	0.27 15.9	3.6	14.8	0.15
		GFPCA method GFPCA method		AGB method AGB method
		Table 5.2: Comparison of pansharpening algorithms for enhanced EDS abundance
		images of dataset 2.

Table 5 .

 5 4: Comparison of pansharpening algorithms for enhanced EDS abundance images on simulated SEM image, for different forcing levels.

	5.7. CONCLUSION

  Let us fix a point x and use the fact that the data points on the training set are i.i.d. that follow P (x), {x i |i ∈ [1, n]} i.i.d. ∼ P (x), to obtain:

	CHAPTER 6. SPATIAL REGRESSION FOR IMAGE PANSHARPENING:
	APPLICATION TO MULTIMODAL SEM IMAGE FUSION	145

σ 1 ) 2 P (x)dx.

  x 1x 2 ) . . . γ θ * (x 1x n ) 1 γ θ * (x 1x 2 )

				
	. . . γ θ * (x 1 -x n ) 1	0 . . . . . . . . .	. . . . . . γ θ * (x n-1 -x n ) . . .        w 1 (x * ) w 2 (x * ) . . . w n (x * )       	γ θ * (x 2 -x n ) 1 . . . . . . 0 1       1 0 =        γ θ  •      	(6.63)
			-µ

* (x 1x * ) γ θ * (x 2x * ) . . . γ θ * (x nx * ) 1

  CHAPTER 6. SPATIAL REGRESSION FOR IMAGE PANSHARPENING: APPLICATION TO MULTIMODAL SEM IMAGEFUSION 155 where const represents a constant that does not depend on ω.By taking now p → ∞, it leads to the following true risk function:E HS|X=x * (w) = var hs j (x * ) -

			n
			w i (x * )hs(x i )	(6.66)
			i=1
			n
		-2λ	j=1
			) =
	const -	σ 2 2

p j=1 (hs j (x * ) -n i=1 w i (x * )hs(x i )) 2 λ n i=1 C r (x * , x i )w i (x * ),

Table 6 .

 6 2: Comparison of spatial interpolation algorithms for enhanced EDS abundance images of dataset 2.

		CC SAM RMSE ERGAS CC g
	AGB	0.30 15.7	1.4	20.1	0.45
	regularized kriging λ = 0.3 0.28 15.7	1.31	49.2	0.23
	regularized kriging λ = 0.6 0.21 60.6	39.8	73.1	0.13

  Nous avons proposé une nouvelle approche pour la classification de pixels d'images multi/hyper-spectrales. Le but de cette technique est de représenter et de décrire de façon efficace les caractéristiques spatiales / spectrales de ces images. Ces descripteurs multi-échelle profond visent à représenter le contenu de l'image tout en tenant compte des invariances liées à la texture et à ses transformations géométriques.

	Cette thèse porte sur la statistique spatiale
	multivariée et l'apprentissage appliqués aux
	images hyperspectrales et multimodales. Les
	thèmes suivants sont abordés :
	Fusion d'images :
	Le microscope électronique à balayage
	(MEB) permet d'acquérir des images à partir
	d'un échantillon donné en utilisant différentes
	modalités. Le but de ces études est
	d'analyser l'intérêt de la fusion de
	l'information pour améliorer les images
	acquises par MEB. Nous avons mis en oeuvre
	différentes		techniques	de	fusion	de
	l'information	des images,	basées	en
	particulier sur la théorie de la régression
	spatiale. Ces solutions ont été testées sur
	quelques jeux de données réelles et
	simulées.
	Classification spatiale des pixels d'images
	multivariées :
	Réduction spatiale de dimensionnalité :
	Nous proposons une technique pour extraire
	l'espace des fonctions en utilisant l'analyse
	en composante morphologiques. Ainsi, pour
	ajouter	de	l'information	spatiale	et
	structurelle, nous avons utilisé les opérateurs
	de morphologie mathématique.

1.4. STATE-OF-THE-ART ON CLASSIFICATION

2.2. BASICS ON MORPHOLOGICAL IMAGE REPRESENTATION(a) (b)

4.2. NOISE ON SEM IMAGES

Traitement de l'image, Machine Learning, Méthodes à noyaux, Morphologie mathématique, Analyse en Composantes Principales, Support Vector Machine, Apprentissage profond, Transformée de scattering, Krigeage.

Remerciements

SUGGESTIONS FOR FUTURE WORK

Bibliography

Chapter 3

Invariant Spatial Classification of Multi/Hyper-spectral Images

Abstract

In this chapter, a novel approach for pixel classification in multi/hyper-spectral images is proposed, leveraging on both the spatial and spectral information of the data. The introduced method lies on a recently proposed framework for learning on distributions of texture descriptors -by representing them with mean elements in reproducing kernel Hilbert spaces (RKHS). These descriptors aim at representing the content of the image while considering invariances related to the texture and to its geometric transformations, so-called spatial invariances. Moreover, we also consider spectral invariances which are related to the physical composition of the pixels. The descriptors are based on the scattering transform, which provides a useful framework for deep learning classification. Moreover, a classification algorithm learning these texture distributions is formulated and interpreted. We also study the consistency and the convergence rate of our classification technique. The performance of the theory for pixels classification is assessed on two hyperspectral images and a set of multispectral images, conventionally used. The results are good in comparison with state of the art.

Résumé

Dans ce chapitre, une nouvelle approche pour la classification des pixels d'images multi / hyper-spectral en s'appuyant à la fois sur l'information spatiale et spectrale des données. La méthode introduite repose sur un ensemble de technique récemment proposé pour pour apprendre la distribution des descripteurs de texture. Ces descripteurs visent à représenter le contenu de l'image tout en considérant les invariances liées à la texture et à ses transformations géométriques. En outre, nous considérons des invariances spectrales qui sont liées à la composition physique des pixels.

Chapter 5

Enhanced EDX images by fusion of multimodal SEM images using pansharpening techniques

Abstract

The goal of this chapter is to explore the potential interest of image fusion in the context of multimodal scanning electron microscope (SEM) imaging. In particular, we aim at merging the backscattered electron images that usually have a high spatial resolution but do not provide enough discriminative information to physically classify the nature of the sample, with energy-dispersive X-ray spectroscopy (EDX) images that have discriminative information but a lower spatial resolution. The produced images are named enhanced EDX. To achieve this goal, we have compared the results obtained with classical pansharpening techniques for image fusion with original approaches tailored for multimodal SEM fusion of information. Quantitative assessment is obtained by means of two SEM images and a simulated dataset produced by a software based on PENELOPE.

Résumé

Le microscope électronique à balayage (MEB) permet d'acquérir des images à partir d'un échantillon donné en utilisant différentes modalités. Le but de ce chapitre est d'analyser l'intérêt de la fusion de l'information pour améliorer les images acquises par MEB. Nous avons mis en oeuvre différentes techniques de fusion de l'information des images, basées en particulier sur des operateurs de mophologie mathematique ainsi que le filtre bilateral. Ces solutions ont été testées sur quelques jeux de données réelles et simulées par un logiciel basé sur PENELOPE.

Introduction

Scanning electron microscope (SEM) is a versatile imaging tool that allows to acquire images with various detectors. Images formed by secondary electrons (SE) reveals mainly topographic contrasts, images formed by backscattered electrons (BSE) in-INFORMATION However, the observed data values hs are corrupted by a Gaussian noise, it is more precise to write:

where is a Gaussian noise of zero mean and standard deviation σ. It is important to notice that the noise is independent of the spatial position. Thus we can rewrite it:

For the sake of simplicity, let us write hs(N 2 ) = (hs(x 1 ), . . . , hs(x n )) .

Moreover, let us consider that we have different observations of the random variable HS at position x * . Finally, let us make the assumption that the data {hs 1 (x * ), . . . , hs p (x * )} are drawn i.i.d. from the previous distribution. Then we have:

where we have considered that ω = (ω 1 (x * ), . . . , ω n (x * )) t . Taking the logarithm of the likelihood function we get:

Therefore, maximizing the likelihood to determine w corresponds to minimizing:

which represents the empirical risk function, which is an empirical variance. In fact, we would like to have the true risk which would corresponds to the limit case when p → ∞. We consider that a good solution is just obtained for all potential data, thus leading us to the following true risk function:

which in matricial form, and using (6.49), it leads to the following problem:

This system of equations is the same than the one of the Gaussian process except for the term µ. This last term is quite important in our case, since it guarantees that the data are on the simplex, which is the case with our abundance maps, and consequently is more physically plausible. But the fundamental difference is that in Gaussian process, the estimated parameters are those optimize the posterior of the regression of the validation set. In geostatistics, the property of stationarity of the field is the key in the optimization.

Let us define a function called the variogram:

which summarizes the variations of the variance of a spatial field according to the translations, i.e., given a distance and a direction. It satisfies the following properties

This new formulation 6.64 does not change the equation system, we just focus on the data near the positions x * , where the locality is added thanks to the backscattered image. We note other works that overcome the issue of using kriging with huge quantity of data [START_REF] Meier | Incremental local gaussian regression[END_REF][START_REF] Nguyen-Tuong | Model learning with local gaussian process regression[END_REF][START_REF] Snelson | Local and global sparse gaussian process approximations[END_REF].

Abstract

This thesis focuses on multivariate spatial statistics and machine learning applied to hyperspectral and multimodal and images in remote sensing and scanning electron microscopy (SEM). In this thesis the following topics are considered:

Fusion of images: SEM allows us to acquire images from a given sample using different modalities. The purpose of these studies is to analyze the interest of fusion of information to improve the multimodal SEM images acquisition. We have modeled and implemented various techniques of image fusion of information, based in particular on spatial regression theory. They have been assessed on various datasets.

Spatial classification of multivariate image pixels:

We have proposed a novel approach for pixel classification in multi/hyper-spectral images. The aim of this technique is to represent and efficiently describe the spatial/spectral features of multivariate images. These multiscale deep descriptors aim at representing the content of the image while considering invariances related to the texture and to its geometric transformations.

Spatial dimensionality reduction:

We have developed a technique to extract a feature space using morphological principal component analysis. Indeed, in order to take into account the spatial and structural information we used mathematical morphology operators