Y Tarabalka

J Chanussot

J A Benediktsson

Segmentation and Classification of Hyperspectral Images Using Watershed Transformation

Keywords: Hyperspectral images, mathematical morphology, watershed, segmentation, classification

Hyperspectral imaging, which records a detailed spectrum of light for each pixel, provides an invaluable source of information regarding the physical nature of the different materials, leading to the potential of a more accurate classification. However, high dimensionality of hyperspectral data, usually coupled with limited reference data available, limits the performances of supervised classification techniques. The commonly used pixel-wise classification lacks information about spatial structures of the image. In order to increase classification performances, integration of spatial information into the classification process is needed. In this paper, we propose to extend the watershed segmentation algorithm for hyperspectral images, in order to define information about spatial structures. In particular, several approaches to compute a one-band gradient function from hyperspectral images are proposed and investigated. The accuracy of the watershed algorithms is demonstrated by the further incorporation of the segmentation maps into a classifier. A new spectral-spatial classification scheme for hyperspectral images is proposed, based on the pixel-wise Support Vector Machines classification, followed by majority voting within the watershed regions. Experimental segmentation and classification results are presented on two hyperspectral images. It is shown in experiments that when the number of spectral bands increases, the feature extraction and the use of multidimen-

Introduction

The growing availability of hyperspectral images has opened the door to numerous new applications in remote sensing and other areas of image analysis.

Hyperspectral sensors capture more than a hundred spectral bands (data channels) simultaneously. Thus, each pixel in a hyperspectral image is presented as the vector of values corresponding to the wide spectrum of reflected light [START_REF] Chang | Hyperspectral Imaging: Techniques for Spectral Detection and Classification[END_REF] (Figure 1 depicts the structure of a hyperspectral image). For instance, NASA Jet Propulsion Laboratory's Airborne Visible-Infrared Imaging Spectrometer (AVIRIS) system has 224 spectral channels with a spectral resolution of around 10nm, covering the wavelengths from 0.4 to 2.5μm [START_REF] Green | Imaging spectroscopy and the airborne visible infrared imaging spectrometer (AVIRIS)[END_REF]. The fine spectral resolution of the data provides an invaluable source of information regarding the physical nature of the different materials, increasing the capability to distinguish structures and objects in the image scene.

However, such a large number of spectral channels implies the high dimensionality of the data and presents challenges to image analysis and classification.

Most of the commonly used methods designed for the analysis of grey level, color or multispectral images are not appropriate for hyperspectral images. As a matter of fact, very limited reference data are usually available (the ratio of the number of referenced samples to the number of spectral channels quickly drops), which limits the performances of supervised classification techniques. Furthermore, for analysis of hyperspectral images the well known curse of dimensionality prevents robust statistical estimations, usual vector norms become meaningless and so on (e.g., the Hughes phenomenon [START_REF] Hughes | On the mean accuracy of statistical pattern recognizers[END_REF]). Therefore, to take full advantage of the rich information provided by the spectral dimension, the development of new algorithms is required.

The first attempts to classify hyperspectral images were designed to assign each pixel to one of the classes based on its spectrum only [START_REF] Landgrebe | Hyperspectral image data analysis[END_REF]. These pixel-level processing systems use a wide range of features, such as the direct spectral information, texture features, and linear and nonlinear transformations of these features. The applied feature extraction procedure often aims at reducing the dimensionality of the data. The features are used for image classification with a wide range of techniques, such as maximum-likelihood or Bayesian estimation techniques [START_REF] Landgrebe | Signal Theory Methods in Multispectral Remote Sensing[END_REF][START_REF] Chen | Statistical pattern recognition in remote sensing[END_REF][START_REF] Bali | Bayesian approach with Hidden Markov modeling and mean field approximation for hyperspectral data analysis[END_REF], neural networks [START_REF] Subramanian | Methodology for hyperspectral image classification using novel neural network[END_REF][START_REF] Yang | A back-propagation neural network for mineralogical mapping from AVIRIS data[END_REF][START_REF] Hernández-Espinosa | Some experiments with ensembles of neural networks for classification of hyperspectral images[END_REF], decision trees [START_REF] Goel | Classification of hyperspectral data by decision trees and artificial neural networks to identify weed stress and nitrogen status of corn[END_REF][START_REF] Zhou | Classification of coastal areas by airborne hyperspectral image[END_REF], genetic algorithms [START_REF] Vaiphasa | Innovative genetic algorithm for hyperspectral image classification[END_REF] and kernel-based methods [START_REF] Camps-Valls | Kernel-based methods for hyperspectral image classification[END_REF][START_REF] Fauvel | Spectral and spatial methods for the classification of urban remote sensing data[END_REF][START_REF] Fauvel | Evaluation of kernels for multiclass classification of hyperspectral remote sensing data[END_REF][START_REF] Guo | Customizing kernel functions for SVM-based hyperspectral image classification[END_REF]. In particular, Support

Vector Machines (SVM) and other kernel-based methods have recently shown good classification results, because they tend to be robust when a limited number of training samples is available.

Further modification to improve classification results consists in the integration of spatial and spectral information in the image analysis. It means that the decision to assign a pixel to a specific class is simultaneously based on the feature vector of this pixel and on some information derived from the pixel's neighborhood. In previous studies, the spatial context was included in the classification system by using morphological filters [START_REF] Fauvel | Spectral and spatial methods for the classification of urban remote sensing data[END_REF], morphological leveling [START_REF] Pesaresi | A new approach for the morphological segmentation of high-resolution satellite imagery[END_REF] and Markov random fields [START_REF] Farag | A unified framework for map estimation in remote sensing image segmentation[END_REF]. These methods show promising results in terms of classification accuracies by incorporating spatial and spectral information. However, these algorithms use fixed-window-based neighborhoods.

This involves the problem of scale selection, especially if the image contains some small or complex structures.

Another approach to define spatial structures consists in performing image segmentation [START_REF] Fu | A survey on image segmentation[END_REF][START_REF] Tarabalka | Segmentation and classification of hyperspectral data using watershed[END_REF][START_REF] Tarabalka | Classification of hyperspectral data using Support Vector Machines and adaptive neighborhoods[END_REF]. The regions obtained from the obtained segmentation map define the spatial context of the pixels within these regions. To make this approach effective, an accurate segmentation of the image is needed.

In previous studies, several methods for multispectral image segmentation have been investigated. Numerous works exploit region merging techniques, where neighboring image segments are merged iteratively based mostly on their spectral similarity. For instance, the eCognition software performs multiresolution segmentation, based on bottom-up region merging [START_REF] Darwish | Image segmentation for the purpose of object-based classification[END_REF]. Initially, each pixel is considered as a separate region, and subsequently pairs of regions are merged, based on a homogeneity criterion, which is a combination of spectral and shape properties. Tilton developed a hierarchical segmentation algorithm [START_REF] Tilton | Analysis of hierarchically related image segmentations[END_REF], which performs region growing and spectral clustering alternately. The main drawback of applying region merging for image segmentation is that the homogeneity criterion, or thresholds must be chosen. For accurate segmentation, these techniques usually produce a pyramid of segmentation maps, using a range of thresholds.

Then, manual interpretation of the results is needed.

Other studies exploit mathematical morphology based segmentation approaches [START_REF] Pesaresi | A new approach for the morphological segmentation of high-resolution satellite imagery[END_REF][START_REF] Soille | Morphological partitioning of multispectral images[END_REF][START_REF] Flouzat | Spatial and spectral segmentation of satellite remote sensing imagery using processing graphs by mathematical morphology[END_REF][START_REF] Scheunders | Multivalued image segmentation based on first fundamental form[END_REF][START_REF] Li | Evaluation of multiscale morphological segmentation of multispectral imagery for land cover classification[END_REF][START_REF] Noyel | Morphological segmentation of hyperspectral images[END_REF][START_REF] Akçay | Automatic detection of geospatial objects using multiple hierarchical segmentations[END_REF], which mostly use granulometries or watershed transformation. The extension of morphological operators to the case of multispectral images is not straightforward, because there is no natural way for total ordering of multivariate pixels, which is a requirement in mathematical morphology. An extensive literature on mathematical morphology for color and multispectral images is available [START_REF] Chanussot | Bit mixing paradigm for multivalued morphological filters[END_REF][START_REF] Chanussot | Total ordering based on space filling curves for multivalued morphology[END_REF][START_REF] Lambert | Extending mathematical morphology to color image processing[END_REF][START_REF] Hanbury | Morphological operators on the unit circle[END_REF][START_REF] Angulo | Morphological coding of color images by vector connected filters[END_REF][START_REF] Aptoula | A comparative study on multivariate mathematical morphology[END_REF]. In particular, the watershed segmentation of color images was investigated in [START_REF] Shafarenko | Automatic watershed segmentation of randomly textured color images[END_REF][START_REF] Chanussot | Watershed approaches for color image segmentation[END_REF][START_REF] Angulo | Mathematical morphology in color spaces applied to the analysis of cartographic images[END_REF].

However, the above morphological and watershed methods are not suitable for segmentation of hyperspectral images, due to the following reasons:

• A hyperspectral image is composed of hundreds of spectral channels.

Therefore, the use of total ordering schemes for multivariate data, such as the bit mixing paradigm [START_REF] Chanussot | Bit mixing paradigm for multivalued morphological filters[END_REF], is not possible, because it would lead to a huge number of values stored for each pixel.

• In previous studies, polar-based representations (HLS, HSV) and perceptional color spaces (LUV, LAB) were used for morphological analysis of color images [START_REF] Angulo | Unified morphological color processing framework in a lum/sat/hue representation[END_REF][START_REF] Shafarenko | Automatic watershed segmentation of randomly textured color images[END_REF]. These methods are not applicable for hyperspectral images.

In a recent paper, Noyel et al. [START_REF] Noyel | Morphological segmentation of hyperspectral images[END_REF] gave an overview of the literature on the watershed-based multispectral image segmentation and performed watershed segmentation of hyperspectral images. Their method is composed of spectral classification to obtain markers and computation of a multivariate gradient to get spatial information. Only visual results (the obtained segmentation maps) are presented in the article. Therefore, the question of defining a watershed transformation for the case of hyperspectral images has only recently been raised in the literature and needs further investigations.

Some studies have been conducted on spectral-spatial classification of multispectral images. Linden et al. [START_REF] Der Linden | Classifying segmented hyperspectral data from a heterogeneous urban environment using support vector machines[END_REF] used the vector mean as a feature for each region in a spectral-spatial classifier. First, they performed the segmentation of a hyperspectral image, based on region growing (using the eCognition software [START_REF] Darwish | Image segmentation for the purpose of object-based classification[END_REF]). Then, a vector mean was computed for each region, such that the value in each band represented the average spectral information of the pixels in this region in the respective band. Afterwards, the regions were classified by an SVM classifier. However, the obtained results were mostly not an improvement over those obtained by the pixel-wise SVM classification. Li and Xiao [START_REF] Li | Evaluation of multiscale morphological segmentation of multispectral imagery for land cover classification[END_REF] used spectral and spatial information for classification of a multispectral (4-band SPOT 5) image. There, a watershed segmentation and a pixel-wise maximum likelihood classification of an image were independently performed.

Then, pixels of the whole region were assigned to one class if more than 50% of pixels in this region were categorized into one class by a pixel-wise classifier.

The classification results were substantially improved with the spectral-spatial approach compared to the pixel-wise maximum likelihood classification.

Widayati et al. [START_REF] Widayati | Application of combined pixel-based and spatial-based approaches for improved mixed vegetation classification using IKONOS[END_REF] used spatial information in a classifier, in order to perform a spectral-spatial classification of a multispectral (4-band IKONOS) image.

First, a segmentation map was obtained using the Merge Using Moments algorithm [START_REF] Cook | MUM (Merge Using Moments) segmentation for SAR images[END_REF]. Then two options were explored. First, each region from the segmentation map was classified using its vector mean as a feature. In another approach, the pixel-wise classification map obtained by maximum likelihood classifier was combined with the segmentation map using majority voting: for every region from the segmentation map, all the pixels were assigned to the majority class within this region. Results of the pixel-wise maximum likelihood classification were also used for comparison. Of all these three methods, the spectral-spatial classification using majority voting gave the highest overall accuracy.

The main contributions of this paper are two-fold:

1. The first contribution is the extension of the watershed segmentation algorithm to hyperspectral data. Watershed transformation is usually applied to the gradient, which must be a scalar function, in order to form a complete lattice structure [START_REF] Soille | Morphological Image Analysis[END_REF]. When processing a hyperspectral image, information about spatial structures must be extracted from all bands in the optimal way. Different approaches to segment the hyperspectral image by watershed are proposed and investigated. In particular, different ways to compute a one-band gradient function of a hyperspectral image are considered. We emphasize that this study is not a comparison of watershed segmentation methods in general. Instead, it focuses on the extension of watershed transformation for hyperspectral images, for the purpose of accurate segmentation and further classification. The paper is organized as follows. In Section 2, an overview of the watershed technique is given, and then the extension of the watershed algorithm to hyperspectral images is discussed. In Section 3, the developed segmentation and classification scheme is presented. In Section 4, experimental results and comparisons are presented and discussed. Finally, conclusions are drawn in Section 5.

Watershed segmentation

Watershed transformation is a powerful mathematical morphology technique for image segmentation [START_REF] Meyer | Morphological segmentation[END_REF][START_REF] Soille | Morphological Image Analysis[END_REF]. It was introduced in image analysis by Beucher

and Lantuéjoul [START_REF] Beucher | Use of watersheds in contour detection[END_REF], and subsequently a lot of algorithms for its implementation have been proposed.

The watershed transform considers a two-dimensional one-band image as a topographic relief. The value of a pixel h stands for its elevation. The watershed lines divide the image into catchment basins, so that each basin is associated with one minimum in the image (see Fig. 2). The watershed transformation A wealth of literature describes techniques for computing the watershed transformation (see for instance [START_REF] Vincent | Watersheds in digital spaces: an efficient algorithm based on immersion simulations[END_REF][START_REF] Meyer | Topographic distance and watershed lines[END_REF][START_REF] Stoev | Rafsi -a fast watershed algorithm based on rainfalling simulation[END_REF][START_REF] Osma-Ruiz | An improved watershed algorithm based on efficient computation of shortest paths[END_REF]). A review of watershed algorithms can be found in [START_REF] Vincent | Watersheds in digital spaces: an efficient algorithm based on immersion simulations[END_REF][START_REF] Stoev | Rafsi -a fast watershed algorithm based on rainfalling simulation[END_REF]. Vincent and Soille [START_REF] Vincent | Watersheds in digital spaces: an efficient algorithm based on immersion simulations[END_REF] have proposed an efficient watershed algorithm using flooding simulations, which has become one of the classical algorithms to compute watersheds.

The output of the watershed transform is a partition of the image composed of regions (sets of pixels connected to the same local minimum) and of watershed pixels (WHEDs, the borders between the regions). Figure 3 shows an example of watershed transformation in one dimension, where three regions, associated with the three minima, are defined. The two maxima correspond to the borders between regions and are not assigned to any region (watershed pixels).

Typically, the result of watershed segmentation on the gradient image without any additional processing is a severe oversegmentation (every single local minimum of the gradient leads to one region). Common ways to reduce the number of local minima are to filter the original image or the gradient function (e.g., area filtering) or to use markers [START_REF] Meyer | Morphological segmentation[END_REF]. The oversegmentation effect can be also corrected using some post-processing, such as merging of similar neighboring regions.

In the next subsection, different approaches for segmentation of hyperspectral images by watershed are discussed and extended.

Watershed segmentation of hyperspectral images

As previously mentioned, the watershed transformation requires as input a one-band image and gives as a result a one-band segmentation map (where each pixel contains the label of the catchment basin or the watershed pixel label).

In this paper, we aim to apply this transformation to a B-band hyperspectral image. Let us consider this image as a set of n pixel vectors X = {x j ∈ R B , j = 1, 2, ..., n} (each pixel is characterized by its spatial location and a vector of spectral values; see Fig. 1). We denote also the image of every spectral band as X λ , λ = 1, 2, ..., B. Different strategies are possible to compute watersheds.

They are summarized in Figure 4.

Before computing a gradient, feature extraction on the original image can be performed, applying one of the transformations such as the Principal Component Analysis (PCA) [START_REF] Richards | Remote Sensing Digital Image Analysis: An Introduction[END_REF][START_REF] Gupta | Wavelet Principal Component Analysis and its application to hyperspectral images[END_REF], Maximum Noise Fraction (MNF) [START_REF] Green | A transformation for ordering multispectral data in terms of image quality with implications for noise removal[END_REF] and Independent Component Analysis (ICA) [START_REF] Wang | Independent component analysis-based dimensionality reduction with applications in hyperspectral image analysis[END_REF][START_REF] Moussaoui | On the decomposition of Mars hyperspectral data by ICA and Bayesian positive source separation[END_REF]. The aim of this step is to obtain either a one-band image or a multi-band image which would contain enough information to distinguish between spatial structures in the image.

If a one-band image with a good distinguishing capability between structures can be obtained, the algorithm for computing a gradient and watershed :

ρ E (Y ) = δ E (Y ) -ε E (Y ). (1)
If, however, at the input of the gradient step we still have a multi-band image, we can proceed in different ways, that can be grouped into the following three categories:

• to compute a vectorial gradient;

• to compute a multidimensional gradient;

• to combine watershed segmentation maps a posteriori.

These three options are discussed in the next three subsections.

Computation of a vectorial gradient

Vectorial gradients are based on the distance between pixel vectors, and produce from the B-band image one-band gradient [START_REF] Evans | A morphological gradient approach to color edge detection[END_REF][START_REF] Noyel | Morphological segmentation of hyperspectral images[END_REF]. Several types of vectorial gradients have been proposed. Noyel et al. [START_REF] Noyel | Morphological segmentation of hyperspectral images[END_REF] proposed to use a metric-based gradient for hyperspectral images, which is described as follows:

For each pixel vector x p , let ψ = [x 1 p , x 2 p , ..., x e p ] be a set of e vectors in the neighborhood of x p (set ψ does not contain x p ). For instance, a four-or an eight-neighborhood (e = 4 or e = 8, respectively) can be used. The metric-based gradient is defined as a difference between the supremum and the infimum of the defined distances between x p and vectors from the set ψ:

∇ MB ψ,d (x p ) = sup i∈ψ {d(x p , x i p )} -inf j∈ψ {d(x p , x j p )}. (2)
Various distances can be used to compute gradient from ( 2) such as Euclidean, Mahalabobis, chi-squared distances [START_REF] Noyel | Morphological segmentation of hyperspectral images[END_REF].

Another type of vectorial gradient is the Robust Color Morphological Gradient (RCMG). This gradient was developed for color images by Evans and

Liu [START_REF] Evans | A morphological gradient approach to color edge detection[END_REF]. Here we investigate the use of the RCMG for hyperspectral images.

For each pixel vector x p , let χ = [x 1 p , x 2 p , ..., x e p ] be a set of e vectors within a structuring element E, which defines the neighborhood of the vector x p , and the set χ contains x p . The Color Morphological Gradient (CMG) is computed as

∇ CM χ,d (x p ) = max i,j∈χ {d(x i p , x j p )}, (3)
i.e., the maximum of the distances between all pairs of vectors in the set χ.

Here, various distances can be chosen. If the Euclidean distance is used, (3) can be rewritten as

∇ CM χ,d (x p ) = max i,j∈χ { x i p -x j p 2 }. ( 4
)
One of the drawbacks of the CMG is that it is very sensitive to noise. To overcome the problem of outliers, the authors of [START_REF] Evans | A morphological gradient approach to color edge detection[END_REF] have proposed to use the RCMG. The scheme to make the CMG robust consists in removing the two pixels that are the furthest apart and then finding the CMG of the remaining pixels. This process can be repeated several times until a good estimate of the gradient is obtained. Thus, the RCMG, using the Euclidean distance, can be defined as

∇ RCM χ,d (x p ) = max i,j∈[χ-REMr ] { x i p -x j p 2 }, ( 5
)
where REM r is the set of the r vector pairs removed. The appropriate value of the parameter r in (5) depends on the chosen structuring element E and the amount of noise present in the image, as discussed in [START_REF] Evans | A morphological gradient approach to color edge detection[END_REF]. When a oneband vectorial gradient is computed, it is used as the input of the watershed algorithm.

Multidimensional gradient methods

Another approach to compute a one-band gradient from the multi-band image consists in considering the B-band image as a set of B one-band images. In this case, the gradients of every spectral band can be computed, using for instance a morphological gradient (see [START_REF] Chang | Hyperspectral Imaging: Techniques for Spectral Detection and Classification[END_REF]). Then the obtained B gradient images ρ E (X λ ), λ = 1, 2, ..., B are combined into one image using linear or non-linear operators.

As an example of the linear operators, the weighted sum of gradients can be computed by

∇ + E (X) = B λ=1 ω λ ρ E (X λ ), (6)
where ω λ denotes the weight of the gradient of the band λ. If ω λ = 1, λ = 1, 2, ..., B, all the bands are supposed to have an equal importance in defining the gradient. Modifying the weight coefficients, the gradient estimation can be improved. For instance, coefficients that are inversely proportional to the estimated noise of each spectral band can be used as the weights in [START_REF] Chen | Statistical pattern recognition in remote sensing[END_REF].

Examples of non-linear operators are the supremum and the median operators. The gradient supremum over morphological gradients of every band is defined as follows: for every pixel the supremum over all gradient images is taken as the output value of this pixel.

Combination of watershed segmentation maps

First, B standard gradients are computed, one for each spectral band. Then each gradient image is used to compute a watershed transformation. This gives B segmentation maps that can be further combined to provide a single output segmentation map.

One of the ways to combine the B segmentation maps, in order to define relevant edges, consists in summing the watershed lines. Here, for each segmentation map obtained from the gradient of band λ, a binary image W λ of watershed lines is produced. Thus, W λ is an image in which watershed pixels are equal to 1 and all other pixels are equal to 0. The sum of watershed lines is computed by

W = B λ=1 W λ . ( 7
)
The obtained image W can be further thresholded, in order to define the border pixels that were presented in most of the segmentation maps, hence ensuring a reliable edge detection.

However, when summing the watershed lines, we do not have information about regions anymore, but only about edges. Furthermore, some edges can become open after thresholding. Therefore, closing of edges and image region labeling must be performed after the procedure described above.

Segmentation and classification scheme

As previously mentioned, the information about spatial structures defined by the watershed segmentation algorithm can be used to improve the results of classification of a hyperspectral image. In this section, a new combined spectral-spatial classification scheme is presented for hyperspectral images based on watershed segmentation.

The general flow-chart of the proposed segmentation and classification scheme is given in Figure 5. At the input we have a B-band hyperspectral image X = {x j ∈ R B , j = 1, 2, ..., n}. First, the watershed segmentation is performed on this image, using one of the approaches discussed in Section 2 and summarized in Figure 4. In the resulting segmentation map, each pixel contains the label of the region it belongs to, or the watershed pixel label (see Figures 3,6).

It is often desirable to produce a segmented image where each pixel belongs to some region, without border pixels between regions. In this case, each watershed pixel can be assigned to one of the regions in its neighborhood. For this pur- pose, we propose to compute for every region S (S = {s j ∈ R B , j = 1, 2, ..., l}, S ⊆ X, with l equal to the number of pixels in the region) the standard vector median [START_REF] Astola | Vector median filters[END_REF]. A standard vector median s V M for a set of pixel vectors is a vector, which fulfills the condition that the sum of the distances between this vector and all the other vectors in the set is minimal (for instance, L 1 norm is used to compute distances):

s V M = arg min s∈S ⎧ ⎨ ⎩ l j=1 s -s j 1 ⎫ ⎬ ⎭ . ( 8
)
Every watershed pixel is assigned to the neighboring region with the "closest" median, i.e., the distance between the vector median of this region and the watershed pixel vector is minimal (see the example in Fig. 6).

After the image is segmented into regions, this spatial information should be used to improve the classification results. Two approaches to integrate spectral and spatial information into the classification system can be distinguished:

1. To define a feature or a set of features for each region from the segmentation map, and classify regions using these features.

2. To perform a pixel-wise classification first, and then combine a pixel-based classified image with the segmentation results.

Here we propose to use the second approach for the spectral-spatial classification of hyperspectral images, in order to improve the results of the pixel-

No WHEDs

1 1 1 1 1 1 1 0 0 0 0 0 1 0 2 2 2 2 1 0 0 2 2 2 1 1 0 0 0 0 0 0 0 3 3 3 3 3 3 3 3 3 1 1 1 1 1 1 1 1 1 2 2 2 1 1 2 2 2 2 1 1 2 2 2

Majority vote within watershed regions

With WHEDs

Results of spectral-spatial classification

No WHEDs With WHEDs

Watershed pixels to the neighboring regions wise classifier using the spatial information from the segmentation map (see Figures 5,[START_REF] Chen | Statistical pattern recognition in remote sensing[END_REF]. The proposed approach is outlined in [START_REF] Tarabalka | Segmentation and classification of hyperspectral data using watershed[END_REF]. First, a pixel-wise classification by SVM on the original hyperspectral image is performed. Then, for every watershed region S, all the pixels are assigned to the most frequent class within this region (majority vote approach). Please note that majority vote within fixed neighborhoods is a standard spatial regularization procedure following a pixel-wise classification. Here, we propose to use the results of a segmentation to define an adaptive neighborhood for each pixel.

The watershed pixels can be either left not processed during the majority voting (what we call the No WHEDs approach), or assigned to the regions with the "closest median" before the majority voting is performed considering all the pixels (what we call the With WHEDs approach). Figure 6 shows an example of the combination of the pixel-wise classification map and the segmentation map by majority vote, using both No WHEDs and With WHEDs approaches.

Experimental results and discussion

Segmentation and classification of the

Experimental results

Different approaches of the watershed transformation of hyperspectral images described in Section 2 were tested. Watershed segmentation was performed on the gradient function obtained in four different ways:

1. Band50 : A morphological Beucher gradient (see [START_REF] Chang | Hyperspectral Imaging: Techniques for Spectral Detection and Classification[END_REF]) was computed on one band. Band No. 50 was chosen arbitrarily, but similar results were obtained with other non noisy bands. Here and in all the following morphological operations, a 3 × 3 square structuring element E was used (center of E was in the center of the square).

2. SumBands: First, a morphological gradient for every band was computed.

Then, the sum of gradients was obtained by [START_REF] Chen | Statistical pattern recognition in remote sensing[END_REF], assuming that ω λ = 1, λ = 1, 2, ..., B.

3. Sum4PCA: First, the PCA transformation was performed on the original image (using the ENVI software [59]). The first four principal components contained 99.16% of the total variance in the data. Morphological gradients of the first four PCA components were computed, and then summed together (using ( 6), with equal unitary weights).

RCMG:

The RCMG on the original image was computed by ( 5), with

r = 1.
The obtained four gradient images are shown in Figure 8. The principal borders of objects are defined by all gradients. Thus, most of the spatial structures can be recognized only from the single band No. 50. However, some structures, for instance shadows were not defined by the gradient Band50 (see Fig. Furthermore, the watershed transformation was applied to each of the four obtained gradients, using the algorithm of Vincent and Soille [START_REF] Vincent | Watersheds in digital spaces: an efficient algorithm based on immersion simulations[END_REF] (based on 8neighborhood connectivity). The resulting segmentation maps for watershed on Band50, SumBands, Sum4PCA and RCMG gradients contained 11641, 10558, 10345 and 11802 regions, respectively. Figure 9.(a) shows the segmentation map based on the RCMG, where the main spatial structures can be seen. (The colors of each region correspond to the label of this region, scaled in order to obtain a grey-scale 8-bit image.)

As expected, the obtained watershed results were severely oversegmented.

Objects were represented mostly by several regions. The first aim was to obtain the segmentation map where each region contained pixels belonging to one object, i.e., where there were no undersegmentation errors. As mentioned in Section 2, oversegmentation can be corrected by merging regions. On the other hand, it may be more difficult to cope with the undersegmentation problem.

To assess qualitatively the accuracy of the watershed segmentation, the obtained segmentation maps were incorporated into the combined spectral-spatial classifier. First, for every segmentation map, the vector median of every region was computed, and the watershed pixels were assigned to one of the neighboring regions based on the minimal distance between the watershed pixel and the vector median of the different regions.

The multi-class pairwise (one versus one) SVM classification, with the Gaussian Radial Basis Function (RBF) kernel, of the original hyperspectral image was performed, using the LIBSVM library [START_REF] Chang | LIBSVM -a library for support vector machines[END_REF]. The optimal parameters C and γ were determined by 5-fold cross validation: C = 128, γ = 0.125. Figure 9.(b)

shows the obtained classification map.

After the pixel-wise SVM classification, the majority vote within the watershed regions was performed. Both the No WHEDs and With WHEDs approaches introduced in Section 3 were applied for each of the four segmentation maps. The global classification accuracies for the pixel-wise SVM and combined spectral-spatial classification are presented in Table 1. Here and in the following the accuracies are based on results for test data. The following measures of accuracy were used:

• Overall Accuracy (OA) is the percentage of correctly classified pixels;

• Average Accuracy (AA) is the mean of class-specific accuracies, i.e., the mean of the percentage of correctly classified pixels for each class;

• kappa coefficient (κ) is the percentage of agreement (correctly classified pixels) corrected by the number of agreements that would be expected purely by chance [START_REF] Richards | Remote Sensing Digital Image Analysis: An Introduction[END_REF].

The number of training and test samples along with the class-specific accuracies for test data are given in Table 2. Figure 9.(c) shows the classification map for the spectral-spatial classification, using the watershed map on the RCMG and the With WHEDs approach. In order to compare the obtained results with previous works that used an SVM and spatial information for hyperspectral image classification, we have included in Tables 1 and2 accuracies of mathematical morphology-based classification of the University of Pavia image using an SVM, principal components and extended morphological profiles (EMP); results are taken from Plaza et al. [START_REF] Plaza | Recent advances in techniques for hyperspectral image processing[END_REF], where the same training and test samples were used for classification. This method was recently proposed by Benediktsson et al. [START_REF] Benediktsson | Classification of hyperspectral data from urban areas based on extended morphological profiles[END_REF] and is considered as one of the most advanced methods for spectralspatial classification of a multi-band data. Other results of joint spectral-spatial classification of the considered image can be found in [START_REF] Aksoy | Spatial techniques for image classification[END_REF][START_REF] Fauvel | Spectral and spatial methods for the classification of urban remote sensing data[END_REF][START_REF] Fauvel | Spectral and spatial classification of hyperspectral data using SVMs and morphological profiles[END_REF]. Classification accuracies were improved by the spectral-spatial classification for almost all the classes (see Table 2). For some classes, like alphalt, bitumen and bricks, accuracies were much improved by including spatial information, i.e., the accuracies improved by 5.81 to 9.72 percentage points, mostly because of the noise reduction. Here, different segmentation approaches led to the best classification accuracies for different classes. The RCMG segmentation map improved the classification accuracies for the classes meadows, trees, metal sheets and bitumen. For the classes alphalt and bare soil, the best accuracies were achieved when applying the Sum4PCA segmentation. That confirms the assumption that the Sum4PCA gradient defined well the information about road structures. For the classes gravel and bricks, the best accuracies were achieved using the Band50 and SumBands approaches, respectively.

Concluding discussion

Based on the above, the following conclusions can be drawn about the accuracy of the watershed segmentation:

1. The RCMG leads to the best segmentation and classification results. Applying the vectorial gradient based on the Euclidean distance between pixel vectors to the hyperspectral image gave the best classification accuracies, despite the high-dimensionality of data.

2. The watershed segmentation based on the gradient of one band is the least accurate approach as was expected, since it may be difficult to distinguish some different neighboring structures when using only one particular band.

Still, as the inclusion of the information from the Band50 segmentation map improved the classification accuracies, most of the spatial structures could be retrieved from this single band.

3. The summing of the gradients of the first four PCA components gave slightly better results than the summing of the gradients of all bands, which indicates that the first PCA components contain the most important spatial information. Thus, applying the feature extraction before computing the gradient and watershed leads to better segmentation results.

Thus, the incorporation of spatial information into the classifier using the majority vote approach led to more homogeneous objects in the resulting classification map, when compared to the pixel-wise classification. However, when performing watershed segmentation, it is usually difficult to identify small but significant structures as separate regions. They may be identified as the border pixels and then be assimilated with one of the neighboring regions. Therefore, the classification accuracies of small classes, in our case trees and shadows, are not improved significantly, or reduced. It can be also seen on the left-bottom part of the image that some small shadows regions were assimilated with the regions in their neighborhood.

When we compare the obtained results with the recent results of spectralspatial classification using SVM and EMPs (see Tables 1 and2), the proposed segmentation and classification approach leads to higher global accuracies. Furthermore, accuracies for 4 from 9 classes are improved by our technique.

Segmentation and classification of the

Experimental results

In this experiment, we do not choose one single band to perform a gradient, because there are more classes, with similar spectral responses, and some classes can be confused using only one specific band. Thus, three gradient functions were computed: SumBands, Sum4PCA, RCMG. They were computed in the same way as for the University of Pavia image, as described in the previous subsection. For the Indiana image, the first four PCA components contain 99.23% of the total variance for the data.

As for the previous data set, the watershed transformation was applied to each of the gradient functions. The segmentation maps contain 1215, 1097 and 1277 regions for the SumBands, Sum4PCA and RCMG gradients, respectively. the information from the segmentation maps was used for the spectral-spatial classification, as described below.

First, the segmentation maps for the With WHEDs approach were computed, by assigning the watershed pixels to the neighboring regions that had the closest Previous classification results for the Indiana image can be found in [START_REF] Tadjudin | Classification of high dimensional data with limited training samples[END_REF][START_REF] Plaza | Recent advances in techniques for hyperspectral image processing[END_REF] for comparison. However, the accuracies in the referenced works are not directly compared with those given in this paper, because different training-test sets are used. Furthermore, in order to evaluate the efficiency of the proposed segmentation and classification scheme using watershed, we have included in Tables 3 and4 the spectral-spatial classification results, obtained by applying another segmentation technique followed by the procedure proposed in this paper (SVM classification and majority vote within the regions from a segmentation map).

This time, the Hierarchical Image Segmentation (HSEG) algorithm [START_REF] Tilton | Analysis of hierarchically related image segmentations[END_REF] was used to segment a hyperspectral image. The NASA Goddard's RHSEG software provides an efficient implementation of this algorithm [START_REF] Tilton | HSEG/RHSEG, HSEGViewer and HSEGReader user's manual (version 1.40)[END_REF], which we used for our investigation. At the initialization step each pixel is considered as one region. The algorithm merges iteratively the most similar adjacent and nonadjacent regions. We used a Spectral Angle Mapper between the region mean vectors as the dissimilarity criterion between regions. The relative importance of merging of non-adjacent regions versus region growing (when only adjacent regions are merged) can be tuned. From experimental results, simple region growing leads to the highest classification accuracies. Segmentation maps for several levels of hierarchy were chosen interactively, and the results leading to the best classification are shown. The class-specific accuracies were improved by using the spatial information in the classification for almost all the classes (see Table 4). After the combined spectral-spatial classification, the classification accuracy was reduced for only one class (oats). The classification using the With WHEDs approach led to 0% accuracy for this class. The oats class has only 20 pixels in the reference map, and represents a small and very narrow (of 2 pixels width) rectangular field.

Thus, after majority voting within the watershed regions, pixels of the oats class were assigned to the classes in its neighborhood (mostly to the grass/trees class). This drawback of spectral-spatial classification, i.e., the risk of loosing small spatial structures, was discussed in the previous subsection.

Concluding discussion

In assessing the watershed segmentation approaches we note that the results for the Indiana image are different from those for the University of Pavia image.

For the Indiana image, the linear combination (sum) of gradients of all bands or several principal components leads to better segmentation results than the use of the vectorial gradient. This can be explained by the fact that the Indiana If we analyze classification results obtained by majority vote within the HSEG regions, we can conclude that:

1. The global and most of class-specific accuracies are improved when compared to pixel-wise classification. This proves the efficiency of the proposed spectral-spatial classification scheme.

2. The accuracies are lower than those obtained when using watershed for segmentation (except for 2 classes). This confirms why it is desirable to use watershed segmentation for hyperspectral images.

Conclusions

Hyperspectral imaging provides rich spectral information per pixel, increasing the capability to distinguish physical structures in a scene. However, a large number of spectral channels presents challenges to image classification. While pixel-wise classification techniques process each pixel independently without considering information about spatial structures, further improvement of classification performances can be achieved by the incorporation of spatial information into classifier, especially in areas where structural information is important to distinguish between classes.

In this paper, there are two main contributions:

1. The extension of the watershed segmentation algorithm for hyperspectral images was proposed. In particular, different ways were investigated to obtain a one-band gradient function from a hyperspectral image and the combination of watershed segmentation maps was proposed.

2. A new methodology, a spectral-spatial classification scheme for hyperspectral images was proposed. The new method is based on the pixel-wise SVM classification, followed by majority voting within the watershed regions.

Thus, segmentation defines an adaptive neighborhood for each pixel.

In two experiments it was shown that the combined spectral-spatial classification, based on majority voting within the regions obtained by the watershed segmentation algorithms, led to higher classification accuracies when compared to pixel-wise classification or previously proposed techniques. Furthermore, classification maps with more homogeneous regions were obtained with the proposed approach.

In conclusion, the proposed classification methodology succeeded in taking advantage of spatial the spectral information simultaneously. Furthermore, experimental results have revealed that it is desirable to use the watershed segmentation map in a spectral-spatial classifier. The watershed transformation is a completely unsupervised method since it does not require any input parameters (thresholds). Therefore, it can be incorporated into an automatic classification system.

The proposed spectral-spatial classification scheme is especially suitable in classifying images with large spatial structures. The drawback of the proposed method is that when applying watershed segmentation, small spatial structures are often not identified as separate regions. This leads to the assimilation of these structures with larger neighboring structures when majority voting is performed within the watershed regions.

In our future work, we will attempt to improve the segmentation results. In particular, we envision to explore the use of additional filtering and merging of regions for that purpose.

Figure 1 :

 1 Figure 1: Structure of a hyperspectral image

2 .

 2 The obtained watershed segmentation map is further incorporated into a spectral-spatial classifier, aiming at improving classification accuracies, when compared to pixel-wise classification. Thus, the second contribution of this paper is the development of the new segmentation and classification scheme to analyze hyperspectral data. The proposed classification method combines results of a pixel-wise SVM classification and the segmentation map using majority vote approach[START_REF] Tarabalka | Segmentation and classification of hyperspectral data using watershed[END_REF]. This contribution can be summarized as follows: the segmentation defines an adaptive neighborhood for each pixel which is used for the spatial regularization following a pixel-wise classification. To test the developed segmentation and classification algorithms, two hyperspectral airborne images are used: A 103-band ROSIS (Reflective Optics System Imaging Spectrometer) image of the University of Pavia, Italy, and a 220-band AVIRIS image taken over the Northwestern Indiana's Indian Pine site [45].

Figure 2 :

 2 Figure 2: Topographic representation of a one-band image

Figure 3 :

 3 Figure 3: Example of watershed transformation in one dimension.

Figure 4 :

 4 Figure 4: Flow-chart which shows strategies of applying watershed to a hyperspectral image

Figure 5 :

 5 Figure 5: Flow-chart of the proposed segmentation and classification scheme.

Figure 6 :

 6 Figure 6: Example of spectral-spatial classification.

  University of Pavia image 4.1.1. Dataset The University of Pavia image is of an urban area that was recorded by the ROSIS-03 optical sensor. The image scene is the urban area surrounding the University of Pavia, Italy. The image has spatial dimensions of 610 by 340 pixels, with a spatial resolution of 1.3m per pixel. The number of bands of the ROSIS-03 sensor is 115 with a spectral coverage ranging from 0.43 to 0.86μm. The 12 most noisy channels have been removed, and the experiments are conducted on the 103-band image. Nine classes of interest are considered, with the number of test and training samples detailed for each class in Table 2. Training-test set was provided by Prof. Paolo Gamba together with the hyperspectral image. False color image of the University of Pavia and the reference data are presented in Figure 7.

Figure 7 :

 7 Figure 7: University of Pavia image: (a) Three-band color composite; (b) Reference data: Asphalt, meadows, gravel, trees, metal sheets, bare soil, bitumen, bricks, shadows and nonlabeled pixels.

Fig. 8 .

 8 Fig. 8.(d)) defines borders that are thinner and clearer than those defined by the other three gradients.

Figure 8 :

 8 Figure 8: Gradients of the University of Pavia image: (a) Gradient of band No. 50; (b) Sum of gradients over all bands; (c) Sum of gradients for the first four PCA components; (d) RCMG using the Euclidean distance, r = 1.

Figure 9 :

 9 Figure 9: University of Pavia image: (a) Watershed on the RCMG, using the Euclidean distance, r = 1; (b) Classification map for the SVM classification; (c) Classification map for the spectral-spatial classification (using watershed map on the RCMG, With WHEDs approach).

  Figures 9.(b) and 9.(c)). The spectral-spatial classification reduces significantly the noise in the classification map, resulting in more homogeneous regions in the output map. It can be also seen on the left-bottom part of the image that some

  Indiana image 4.2.1. Dataset In the second case study, the developed segmentation and classification algorithms are tested on a hyperspectral image of a rural area (the Indiana image) with more bands and a lower spatial resolution as compared to the University of Pavia image. The Indiana image was captured by the AVIRIS sensor over the Indian Pine test site in Northwestern Indiana [45]. The image is 145 by 145 pixels, and the spatial resolution is 20m per pixel. It is composed of 220 spectral channels, and the full spectral range was used for experiments. The data contains 16 classes, which are detailed in Table 4, with a number of samples for each class in the available reference data. Figure 10 depicts the 3-band false color composite and the reference data. In our experiments, we have chosen randomly 10% of the samples for each class from the reference data as training samples (in order to test classification performances when limited reference data are available). (a) (b)

Figure 10 :

 10 Figure 10: Indiana image: (a) Three-band color composite (bands 50, 27, 17); (b) Reference data: Corn-no till, Corn-min till, Corn, Soybeans-no till, Soybeans-min till, Soybeans-clean till, Alfalfa, Grass/pasture, Grass/trees, Grass/pasture-mowed, Hay-windrowed, Oats, Wheat, Woods, Bldg-Grass-Tree-Drives, Stone-steel towers and non-labeled pixels.

Figure 11 .Figure 11 :

 1111 Figure 11.(a) shows the segmentation map obtained from the SumBands gradient. It is difficult to evaluate the accuracy of segmentation from the image, since it is strongly oversegmented. In order to perform such an evaluation,

  vector median. Then, the multi-class one versus one SVM classification, with the Gaussian RBF kernel, of the original Indiana image was performed. The parameters C and γ were computed by 5-fold cross validation giving: C = 1024, γ = 2 -7 . The resulting classification map is given in Figure 11.(b). The classification accuracies for test samples are shown in Tables 3 and 4. The low spatial resolution of the Indiana image leads to the presence of highly mixed pixels which complicates the classification problem. Furthermore, some classes represent small crop fields, and the number of samples in the reference data for different classes varies from 20 to 2468 samples per class. To train the SVM classifier, 10% of the samples for each class were chosen randomly from the reference data. Therefore, some classes were represented by only a few samples in the training set (2 samples for the class oats), which may not provide sufficient representatives for the small classes. Despite these problems, 78.76% of the pixels in the test set were correctly classified by the SVM classifier. As expected, the classification accuracies for the classes that were represented by only a few training samples are low (for instance, only 22.22% accuracy was observed for the class oats). The low classification accuracies for the classes alfalfa, grass/pasture-mowed and oats reduced the average classification accuracy to 69.66%. Now we turn to the spectral-spatial classification. As in the first experiment, the watershed segmentation maps were used in the spectral-spatial classification based on the majority vote method, with the No WHEDs and With WHEDs approaches. Tables 3 and 4 give global and class-specific classification accuracies for the spectral-spatial classification of the Indiana image, respectively. Figure 11.(c) shows the classification map for the spectral-spatial classification, using the watershed map on the SumBands gradient and the With WHEDs approach.

  image contains more spectral channels (220 channels versus 103 channels for the University of Pavia image). Furthermore, the vectorial gradient, based on the Euclidean distance between pixel vectors does not give accurate results due to the curse of dimensionality. In this case, it appears to be more appropriate to compute marginal gradients of every band and to sum them together or to reduce the spectral dimension by performing feature extraction first, and then compute a gradient. The spectral-spatial classification improved classification accuracies when compared to pixel-wise classification. The improvement of the global accuracies was more significant than for the University of Pavia image. From a visual observation (see Figures 11.(b) and 11.(c)), the classification map obtained by the spectral-spatial classification is seen to be much less noisy than the one obtained by the pixel-wise classification. From the obtained results, we can draw the conclusion that the Indiana image contains a lot of relatively large regularly structured spatial objects (crop fields) which is the reason why the inclusion of the spatial information was very efficient and improved the classification accuracies significantly.

Table 1 :

 1 Global classification accuracies in percentage for the University of Pavia image: Overall Accuracy (OA), Average Accuracy (AA) and kappa coefficient (κ).

	Method		OA	AA	κ
	Pixel-wise SVM		81.01	88.25	75.86
		Band50	83.00	89.59	78.31
	No WHEDs	SumBands 83.10 Sum4PCA 83.36	89.80 89.85	78.44 78.76
	SVM +	RCMG	83.48	90.03	78.89
	Majority Vote	Band50	84.83	90.68	80.57
	With WHEDs	SumBands 85.04 Sum4PCA 85.35	90.73 91.09	80.83 81.23
		RCMG	85.42 91.31 81.30
	EMP		85.22	90.76	80.86

As can be seen from Table

1

, all the SVM-based classifiers give high classification accuracies. Furthermore, the combining of spatial information obtained by watershed segmentation with the spectral classification results improves substantially the classification accuracies. The With WHEDs approach gave the best accuracies for all four segmentation results. The best global accuracies are achieved when performing the spectral-spatial classification based on the segmentation map on the RCMG and applying the With WHEDs approach. In that case, the overall accuracy is improved by 4.41 and the average accuracy is improved by 3.06 percentage points compared to the pixel-wise SVM classification. The integration of the different segmentation results in the classification improves the accuracies differently: SumBands performs better than Band50, Sum4PCA better than SumBands, and RCMG performs best.

Table 2 :

 2 Information classes, training-test samples, and class-specific accuracies in percentage for the University of Pavia image.

		EMP		95.36	80.33	87.61	98.37	99.48	63.72	98.87	95.41	97.68
			RCMG	93.64	75.09	66.12	98.56	99.91	97.35	96.23	97.92	96.98
		With WHEDs	SumBands Sum4PCA	93.62 94.65	74.49 74.50	64.79 67.05	97.18 97.97	99.64 98.83	97.83 98.14	94.29 93.99	98.78 98.42	95.97 96.23
	SVM + Majority Vote		RCMG Band50	89.42 92.64	73.33 74.45	66.72 68.76	97.94 95.88	99.73 99.46	94.93 97.31	93.88 93.37	96.20 98.42	98.11 95.85
		No WHEDs	SumBands Sum4PCA	89.70 89.61	72.51 72.83	67.55 67.82	97.36 98.15	99.73 99.37	94.62 95.71	93.37 92.46	96.49 96.11	96.86 96.60
			Band50	88.50	72.73	68.21	96.84	99.64	94.71	92.05	96.82	96.86
		SVM		84.93	70.79	67.16	97.77	99.46	92.83	90.42	92.78	98.11
	Samples	Train Test	548 6304	540 18146	392 1815	524 2912	265 1113	532 4572	375 981	514 3364	231 795
		Class		1-asphalt	2-meadows	3-gravel	4-trees	5-metal sheets	6-bare soil	7-bitumen	8-bricks	9-shadows

Table 3 :

 3 Global classification accuracies in percentage for the Indiana image: Overall Accuracy (OA), Average Accuracy (AA) and kappa coefficient (κ).As in the previous experiment, the combined spectral-spatial classification improved the classification accuracies, and the With WHEDs approach gave the best performances. The use of the watershed map based on the SumBands gradient led to the highest overall accuracy (an improvement by 15.02 percentage points compared to the pixel-wise classification) and kappa coefficient, while the highest average accuracy was achieved by using the Sum4PCA gradient (it

		Method		OA	AA	κ
		Pixel-wise SVM		78.76	69.66	75.73
		SumBands 87.03	75.24	85.18
	SVM +	No WHEDs	Sum4PCA 86.92 RCMG 86.79	77.79 79.36	85.05 84.88
	Majority Vote	SumBands 93.78 80.53 92.88
			Sum4PCA 91.67 80.61 90.47
		With WHEDs	RCMG	92.48	77.26	91.39
		HSEG		92.20	80.49	91.06

is 10.95 percentage points higher than for the pixel-wise classification). The watershed based on the RCMG gradient led to lower classification accuracies than the one based on the SumBands gradient. And the RCMG approach gave a higher overall accuracy but a lower average accuracy than the Sum4PCA approach.

Table 4 :

 4 Information classes, number of labeled samples, and class-specific accuracies in percentage for the Indiana image.

		HSEG		85.82	84.82	91.00	87.86	98.69	89.87	2.04	94.42	98.37	95.83	99.32	00.00	98.43	98.45	75.73	87.21
			RCMG	87.45	88.42	81.04	83.26	98.47	94.39	44.90	94.42	99.11	00.00	98.41	00.00	99.48	97.68	84.21	84.88
	SVM + Majority Vote	With WHEDs	RCMG SumBands Sum4PCA	83.89 90.24 87.84	74.83 90.68 83.49	70.14 86.73 87.20	80.53 89.56 84.98	89.83 99.33 97.44	84.45 95.48 94.21	46.94 44.90 44.90	90.63 93.53 89.96	96.43 99.11 98.22	100 45.83 50.00	97.51 99.32 97.05	5.56 00.00 00.00	94.76 99.48 99.48	94.59 97.51 95.19	74.85 75.44 80.99	84.88 81.40 98.84
		No WHEDs	Sum4PCA	80.94	78.16	78.20	80.62	90.33	87.70	46.94	91.52	94.80	37.5	98.19	16.67	95.29	92.70	76.32	98.84
			SumBands	84.04	78.03	71.56	81.19	89.11	88.25	24.49	92.19	97.33	45.83	98.41	5.56	95.81	93.73	74.56	83.72
		SVM		74.59	64.58	58.77	69.76	79.21	75.41	32.65	87.05	92.72	29.17	96.37	22.22	90.58	91.07	65.50	84.88
	No. of samples	in reference data	1434	834	234	968	2468	614	54	497	747	26	489	20	212	1294	380	95
		Class		1-corn-no till	2-corn-min till	3-corn	4-soybeans-no till	5-soybeans-min till	6-soybeans-clean till	7-alfalfa	8-grass/pasture	9-grass/trees	10-grass/pasture-mowed	11-hay-windrowed	12-oats	13-wheat	14-woods	15-bldg-grass-tree-drives	16-stone-steel towers

Acknowledgments

This research is supported in part by the Marie Curie Research Training Network "HYPER-I-NET". The authors would like to thank Paolo Gamba and David Landgrebe for providing the hyperspectral data and James Tilton for providing the NASA Goddard's RHSEG software.


Nazim Fatès

Nikolaos Vlassopoulos

A Robust Aggregation Method for Quasi-blind Robots in an Active Environment ⋆

Keywords: Decentralized Gathering, Agent Aggregation, Swarm Robotics, Reaction-diffusion-chemotaxis

The aggregation of a swarm of autonomous agents into compact clusters is often a required behavior of multi-agent systems. In the case where no central control or coordination exists, this problem is known as the Decentralized Gathering. This presents a first step in the application of a bio-inspired aggregation scheme to aggregate robots whose abilities are very restricted: the only way they can communicate is through an active environment (stigmergy) and the only information they can receive is the local detection of the waves produced by other robots. The active environment obeys a cellular automaton rule and is simulated with a projection of light on the robots.

Introduction

Let us consider a system constituted of a great number of autonomous robots randomly spread on a plane area, the environment. The robots need to gather to realize a given task, for instance to exchange of information about the area or to unite to collectively remove an obstacle. In the case where the agents are not aware of their position, orientation or of the presence other agents in the environment and if furthermore there are no means of centralized coordination, the problem is known as Decentralized Gathering.

Similar problems have been studied before, using different assumptions about the agents and the environment. An early work [START_REF] Sugihara K | Distributed Algorithms for Formation of Geometric Patterns with many Robots[END_REF] on the subject deals with the distributed coordination of a swarm of robots in order to form geometrical shapes. The work presented in [START_REF] Garnier S | The Embodiment of Cockroach Aggregation Behavior in a Group of MicroRobots[END_REF][START_REF] Garnier S | Self-organized aggregation triggers collective decision making in a group of cockroach like robots[END_REF] studies the emergent aggregation of robots that mimic the behavior of Blattela germanica, while [START_REF] Ando | Distributed memoryless point convergence algorithms for mobile robots with limited visibility[END_REF][START_REF] Flocchini | Gathering of Asynchronous Robots with Limited Visibility[END_REF] study the aggregation problem under the assumption of limited visibility. Finally, [START_REF] Prencipe | Impossibility of Gathering by a set of Autonomous Mobile Robots[END_REF] presents some recent theoretical developments on the problem of decentralized gathering.

One approach to solving the decentralized gathering problem consists of imitating the amoeba species Dictyostelium discoideum [START_REF] Fates | Solving the decentralised gathering problem with a reaction-diffusionchemotaxis scheme -Social amoebae as a source of inspiration[END_REF]. The global behavior of this aggregation scheme has been studied in the scope of cellular automata (CA) modeling. However, the scheme is more generic, as it can be considered as a generic behavioral description. The key ideas that characterize the scheme is the existence of an active environment that conveys simple messages between the agents, here called the virtual amoebae. The agents interact with the environment by either initiating the transmission of a message or detecting the existence of messages in their local environment. These two ingredients form the basis for a stigmergic behavior.

Initially, the aggregation scheme involved agents that were modeled either as particles residing on the cells of a CA or as properties of the CA cells [START_REF] Fates | Solving the decentralised gathering problem with a reaction-diffusionchemotaxis scheme -Social amoebae as a source of inspiration[END_REF][START_REF] Spicher | Translating Discrete Multi-Agents Systems into Cellular Automata: Application to Diffusion-Limited Aggregation[END_REF]. The rules that govern the evolution of the cellular automaton mimic reactiondiffusion waves, namely they follow the Greenberg-Hastings CA [START_REF] Greenberg | Pattern formation and periodic structures in systems modeled by reaction-diffusion equations[END_REF].

The first advantage of this scheme is that both agents behavior rules and the active environment are simple and straightforward to implement. Numerical simulations have further shown that the aggregation scheme is robust under different assumptions and perturbations in the environment, i.e. presence of obstacles and noise in the movement of the agents. In this paper, we present our first results on the application of the aggregation scheme to the case of robotic agents and we examine two different series of experiments that use physically different implementations of the same active environment.

The first approach we will present is based on the ALICE micro-robots [START_REF]Autonomous Mini-Robot ALICE[END_REF]. These robots feature a compact design but have a restricted set of abilities: the sensors of the robots we handled were only two sensors that could measure the light intensity of a source located above the robots. The second approach uses the a modified version of the KheperaIII robots, in conjunction with the interactive table ROMEA that has been specifically developed for conducting robotic experiments [START_REF] Charpillet | Interactive table to study interactions between swarms of robots and active environments[END_REF]. These robots have a far more sophisticated design as they interact with the interactive table. Their standard sensors include peripheral infrared transceivers that can be used to measure distance from obstacles and ultrasound distance sensors. Further, they are equipped with color sensors that are capable of measuring the intensity of the (Red, Green, Blue) color components of the surface beneath the robot, which allows them to read the values projected on the interactive table, as well as LEDs that allow the table to track their position. However, as we will describe later, we made a minimalistic use of these features: our goal is indeed to present a simple method for a aggregating the robots, in particular we want the experiment to be set rapidly without having to tune a great range of parameters. This paper is organized as follows: Section 2 presents the aggregation scheme while Section 3 describes the implementation on the two different types of robot agents. Following, Section 4 presents the advantages and drawbacks that we observed during our experiments, as well as the differences we observed between the two types of robots and environment implementations. Finally, we discuss the results in Section 5.

The virtual amoebae aggregation scheme is based on the interaction of virtual agents with an active environment. To implement the active environment where the robots evolve, we use a CA approach in the form of the Greenberg-Hastings (GHCA) model for reaction-diffusion processes. There are two main advantages that make this type of environment attractive: (a) it is easy to implement and simulate in software and, (b) the reaction-diffusion wavefronts travel over arbitrary large distances in the environment without any attenuation. From a practical viewpoint, making the agents aware of the environment and enabling them to interact with it, i.e. linking the agents with the environment, requires a different approach for the two types of robots. We now describe shortly the GHCA and the parameters that govern the environment.

Description of the CA Environment Layer

Let L = {1, . . . , L} × {1, . . . , L} be a two dimensional array of cells where L determines the environment size. Each cell c = (c x , c y ) is associated with a state σ such that σ ∈ {M, . . . , 0} We will call M the excited state, 0 the neutral state and states M -1, . . . , 1 the refractory states. The state of a cell at time t is denoted by σ t c The time evolution of the state of each cell is a function of its current state and the states of the cells in its neighborhood. If we denote by

E t c = {c ′ ∈ N c : σ t c ′ = M }
the number of excited cells in the neighborhood of c at time t, then the state evolution of a cell, c is described by the following set of equations:

σ t c =    M , σ t c = 0 and |E c | > 0 σ t c -1, σ t c ∈ {1, . . . , M -1} 0, otherwise
Our experiments used two types of neighborhood: 4-connected, 8-connected neighborhood and, in the case of the KheperaIII robots, a special circular neighborhood to simulate the diffusion of isotropic waves. Finally, we assumed that the environment boundaries are free, i.e. wavefronts that reach the boundaries are simply absorbed. This is a desired behavior, since, as we will see, we want the agents to move towards the direction of the source of a wavefront.

The reaction-diffusion process can be quantitatively described as waves that start from single excited points and expand outwards. Figure 1 shows the evolution, after t = 20 steps steps of 5 initial wave departures.

Coupling of the Virtual Amoebae with the environment

The basic model assumptions for the behavior of the virtual amoebae (agents) can be summarized in the following:

-Each agent is independent. -There is no direct communication among agents.

-Each agent interacts only with its immediate environment (local interaction). -Agents are not aware of their position, orientation or relative positions of other agents in the environment (quasi-blindness).

With the above in mind, the behavior of each agent can be described as follows:

1. If a reaction-diffusion wavefront is detected in the immediate neighborhood of an agent, and the agent is not on a refractory cell (or in a refractory state) then the agent moves towards the wavefront. 2. If nothing is detected, then the agent initiates a reaction-diffusion wave by modifying the cell of the environment where it resides. This happens with probability p F so that the average number of waves emitted by an agent are λ waves sec . The agent then becomes idle, which means that becomes insensitive to the environment. Therefore, quantitatively, each agent either moves towards a wavefront, if is detects one excited cell in its neighborhood, or, initiates a reaction-diffusion process with a certain probability that we will call the firing rate. The above rules describe the behavior of agents on a discrete lattice [START_REF] Fates | Solving the decentralised gathering problem with a reaction-diffusionchemotaxis scheme -Social amoebae as a source of inspiration[END_REF]. Naturally, their implementation migration of this behavior to agents operating in a physical environment required some changes. The most important problem is to prevent an agent from detecting and following a wave he has just emitted. In the following section we look into the details of this implementation.

The main differences between the original CA-based instance of the aggregation scheme and its robot counterpart concern space and time. In simulations, both the agents and the environment operate in discrete time but in the physical experiments the movements of the robots are continuous. Moreover, in the discrete model, the actual size of the agents is neglected and each agent is simply assigned to one cell as a location. In the physical experiment, the agents have a non-negligible size and they cover many environment cells. As a consequence, the relative length of the environment cells with respect to the agent dimensions, as well the sampling time, that is, the time unit which separates consecutive CA steps, become important, since they affect how the agents perceive and interact with the environment.

Chemotaxic behavior

Before going into further details, we call the reader's attention on the fact that it is the same behavior that is used for the two robots. This behavior was mainly chosen to cope with the limited abilities of the sensors of the ALICE robots. An overview of these robots is shown in Figure 2. The robots are equipped with two wheels, so that controlling their speed and direction can be accomplished by setting the speed of each wheel. The two light sensors are located on the front of the robot. These sensors are connected to analog-to-digital converters, and the intensity of the incident light is measured in a normalized scale from 0 to 100. The wheel motors can be controlled independently by setting their rotation speed. From an agent-based point of view, the behavior of an agent can be described as a control system with two input variables and two output variables: the light sensor measurements and the motor speeds. The main problem to achieve the aggregation is that by using only two sensors, we cannot determine the exact direction of an incoming wave. An example of such ambiguity is shown in Fig. 3: the robot can not distinguish between the two incoming wavefronts. Indeed, from the different time of arrivals of the wavefronts, the agent can only determine if the wave is more or less parallel to the line that connects the two sensors. In other words, in the best case we are able to determine the absolute value of the angle of the incoming waves, but not its sign. One solution to this problem would be to use three or more sensors. Nevertheless, this would require designing a specialized extension module, as well as some extra processing overhead in order to determine both the angle and the sign. The simplest method to overcome the directional ambiguity is, surprisingly, to completely ignore it. As we can observe in Fig. 4 by programming the robots to turn for an angle 90 o towards the direction of the incident wave, i.e. towards the sensor that detects an incoming wave, the robot will eventually turn towards the source of the waves, although this might require more than one waves to reach the robot.

The second major problem is related to the capability of the agents to follow the gradient of the incident reaction-diffusion waves. Imagine that a reactiondiffusion wave, whose intensity gradually varies from a maximum value to zero is detected first from one of the sensors and following from the second one. Assume that the agent eventually turns so that its direction is perpendicular to the incoming wavefronts. The problem that arises is that since there are only two sensors and they are placed perpendicular to the movement direction of the agent, they are both parallel to the wavefronts. As a consequence so they become unable to accurately detect the gradient, since the light intensity varies in the same way on the two sensors. This problem appears both when the robot faces towards or away from the reaction-diffusion waves. To resolve this ambiguity, we used again the simplest possible approach. Instead of trying to follow the gradient of the waves, moving towards the source of the wavefront, we programmed the robots to only react to the presence of wavefronts, that appear as a sudden change in the detected luminosity. When such a change is detected, the robot turns towards the sensor that detected the change. The effect of this approach Fig. 4. Elementary chemotaxis movements is that the agents follow a "zig-zag" track, as it is shown in Fig. 5 towards the source of the waves. Note that, contrary to the ALICE robots, the KheperaIII robots have an adequate number of color sensors [START_REF] Spicher | Translating Discrete Multi-Agents Systems into Cellular Automata: Application to Diffusion-Limited Aggregation[END_REF], so that they can properly detect the exact direction of incoming wavefronts as well follow the gradient of the waves. However, in order to be able to compare the two experimental setups, we used the exact same method in both cases. As we will discuss in the following sections, this approach proved to be robust and efficient despite its simplicity.

Experimental Setup for the ALICE Robots

As we have already seen, ALICE robots are equipped with two light sensors that are located on top of the robot. This arrangement was the basic guideline on the setup of the aggregation experiments, since the reaction-diffusion waves should be projected on the surface where the robots are located from above. We used a setup that was based on a camera and a projector: the images from the camera were processed in real-time in order to identify the positions of the robots on the environment and map them to cell locations in the CA array. By knowing the agent locations and the state of the array, we simulated the environment using FiatLux CA simulator [START_REF] Fates | FiatLux CA Simulator in Java[END_REF] and projected the resulting images back on the surface where the agents were located. Figure 6 shows a snapshot of one experiment, where we can see the agents and the projected waves.

Problems of the experimental setup

The main problem of this setup is the fact that although the agents can "read" their environment and detect incoming wavefronts, it is not possible for them to "write" directly on the environment, i.e. to initiate a wave. In order to overcome this problem, we delegated the wave initiation part of the agent logic to the environment. However, since the fire decisions are not taken by the agents, it becomes impossible for them to know whether they should enter a "refractory" state and ignore their own emitted waves. The solution to this problem was again implemented in the environment, and more precisely by generating a "mask" that prevented from painting the first steps of the evolution of the reaction-diffusion process. Fig. 6. Snapshot of an experiment with ALICE robots, as taken from the camera that detects the robot positions.

Experimental Setup for the KheperaIII Robots

In the case of the KheperaIII robots, the experimental setup was eased, since the interactive ROMEA table and the existing programming interfaces provided all the tools required to track the position of the robots and simulate the reactiondiffusion process. In that case, in order to prevent the robots from detecting the waves they emitted, we introduced a refractory delay after each fire event, during which the robot became insensitive to the presence of waves in the environment. That delay was empirically set as a function of the environment simulation step time, which determines the speed of waves and the positions of the sensors on the robot. Fig. 7. Snapshot of an experiment using the Khepera robots on the interactive ROMEA table. The experiment uses a circular neighborhood with radius 2.5 and M=6

Results and Observations

The most important feature that became apparent as soon as we started running the experiments with the ALICE robots is the fact that the aggregation scheme is rather robust with respect to perturbations in the agent movement. More precisely, although the robots that we used for our first experiments were calibrated so as to turn by 90 o , this was not the case during the experiments and as a result, the robots would turn from 0 to sometimes much over 120 o . However, this did not obstruct the aggregation process and in all experiments we performed, the robots managed to form a single cluster. To study this robustness in more details, we implemented a simplified model in software that modeled as closely as possible the ALICE robots, with configurable sensor and wheel positions as well as turning radius. Using this model we ran several experiments while varying the angle a robot turns each time it detects a wave and also trying out some configurations where the angle was a uniformly distributed random number between 0 and 90 o . In all of the above cases, the robots managed to converge, with the fastest average convergence time obtained for a turning angle of 90 o .

One major difference between the ALICE and KheperaIII experiments was that in the first case we did not implement any method for avoiding obstacles. Although this resulted in the ALICE robots blocking and pushing each other, still, this behavior would only delay the aggregation. On the other hand, with the KheperaIII robots, using obstacle avoidance is mandatory since the robots cannot push each other. This does not hinder the aggregation but sometimes produces deadlock configurations in the experiments where the environment contains static obstacles.

One further advantage of the aggregation scheme is that it can be described with a small number of parameters. More precisely, the main parameters to set the agents are the fire rate and the refractory time, if we consider the turn angle and movement time constants. Regarding the environment, the main parameters are the number of states M and the discretization time, that is, the time between two simulation steps. Regarding the latter, it is sufficient to set it such that the final speed of the reaction-diffusion waves are greater than the turn speed of the robots and greater than the ratio of robot dimensions divided by the refractory time. This allows us to ensure that no wavefronts will be detected more than once.

Two basic features of the aggregation scheme is that it works using only local operations and stigmergy, as well as the fact that it delegates most of the communication overhead to the environment. The operations are local both in terms of the agents but also, in the case of a CA based environment, in terms of the environment. Therefore the aggregation is performed without any need for global communication. The second feature makes the greater part of the system complexity hidden in the environment. This significantly simplifies the agent design and agent logic, but also requires the of active cells and, so far, we do not how to implement this feature physically.

For the robot control, the behaviour is described with a purely discrete state machine, without having to reside to a "control system" approach. This is an advantage since we have a simple design, but on the other hand, it would be difficult to extend this design to more complex environement types. We thus leave for future work the analysis of other types of obstacle avoidance.

Regarding the disadvantages of our method, the presence of static obstacles, for instance blocked agents, seems to create some deadlock configurations. These configurations appear in the cases where an agent is for example in the corner of an obstacle, either concave or convex. So far we have tested several methods for obstacle avoidance, and although we have some good results, the aggregation times are significantly increased.

One other thing to mention is that the agents seem to destabilize in the cases where we have a robot that would emit waves at a high frequency. In this case, the agents can not separate the consecutive waves, and therefore, turn randomly. This situation was prevented here by the use of idles states for robots. However, it could be an important problem in really noisy environments, in the sense of environments detecting many "false" emissions from agents.

The experiments confirm that the aggregation scheme exhibits good properties of avoidance of obstacles. Further, since the scheme is not based on a gradient method, but only on the detection of wavefronts, that is, abrupt changes on the environment, the number of states of the GHCA can be as low as 3 states (M = 2). Our observations indicate that the system's robustness is linked to the fact that it does not depend on the detection of a gradient, but it works more in an "all-or-nothing" fashion. It also appears that the type and size of neighborhood of the GHCA is not so important, as all three neighborhoods that we have tested gave the same or similar results.

A source of improvement of our approach is to improve the aggregation time by minimizing th e number of colliding waves annihilate. Indeed, not all the emitted waves reach all agents and therefore this somehow slows down the aggregation. From a previous research on the aggregation scheme [START_REF] Vlassopoulos | How fast can Virtual Amoebae aggregate? Analysis for the Optimal Firing Rate in an Instance of the Reaction-Diffusion-Chemotaxis Aggregation Scheme[END_REF], we know that there is an optimal value of the firing rate. A future direction of research is thus to see how this optimal value can be determined for the case of the robotic agents.

Conclusions

We have presented some preliminary results on the application of a decentralized aggregation scheme on two different types of agents. This method is based on local interactions of the agents with the environment, and in the case where the environment is modeled through a cellular automaton, the locality is preserved within the environment as well. The main characteristics of the scheme, as we have deduced from our experiments, is that it is robust to perturbations on the movements of agents. The required complexity on the agents is minimal, as well as the required sensor precision. The main difficulty for implementing the aggregation scheme is the design of an active environment. When this environment contains obstacles, one also need to solve some "deadlock" configurations that appear.

Our first results are available through the AMYBIA project webpage1 , where we have included one of the early videos of our ALICE experiments. Several recordings of our KheperaIII experiments are also available on the webpage and are accessible through an applet2 that "replays" the recorded robot positions and orientation. The list of experiments shown there is by no means final and it will expand as we run more experiments to study and evaluate the properties of the scheme and test new methods for solving the problems that we observe.

Fig. 1 .

 1 Fig. 1. Reaction-diffusion process with M=4, t=20 and 5 initially excited cells. Excited cell are marked with bright red color, while cells in the neutral state are marked with black.

Fig. 2 .

 2 Fig. 2. ALICE robot

Fig. 3 .

 3 Fig. 3. An example of the directional ambiguity, where the robot cannot distinguish between the two incident wavefronts

Fig. 5 .

 5 Fig. 5. Zig-zag approach for moving towards the source of reaction-diffusion waves

http://www.loria.fr/˜fates/Amybia/project.html

http://www.loria.fr/ ~fates/Amybia/html/showViewer.html

Acknowledgments

The authors would like to thank Guy Theraulaz of Centre de Recherches sur la Cognition Animale for providing us with the ALICE robots for our experiments.

⋆ This research was funded by an INRIA grant to project AMYBIA