G Noyel
email: noyel@cmm.ensmp.fr

J Angulo

D Jeulin

A new spatio-spectral morphological segmentation for multispectral remote sensing images

A general framework of spatio-spectral segmentation for multispectral images is introduced in this paper. The method is based on classification-driven stochastic watershed by Monte Carlo simulations, and it gives more regular and reliable contours than standard watershed. The present approach is decomposed into several sequential steps. First, a dimensionality reduction stage is performed using Factor Correspondence Analysis method. In this context, a new way to select the factor axes (eigenvectors) according to their spatial information is introduced. Then a spectral classification produces a spectral pre-segmentation of the image. Subsequently, a probability density function (pdf) of contours containing spatial and spectral information is estimated by simulation using a stochastic watershed approach driven by the spectral classification. The pdf of contours is finally segmented by a watershed controlled by markers coming from a regularization of the initial classification.

Introduction

Multispectral (or hyperspectral) images which are composed of several tens or hundreds of spectral bands are nowadays used currently in remote sensing. These spectral bands bring a lot of information concerning the structure of the ground, the vegetation, buildings, etc. However, in order to efficiently process these large amount of data, one has to develop new methods to analyse these images and especially to segment them (i.e. to group similar pixels into connected classes). Standard methods of segmentation require to be extended to this type of rich images.

Two types of information are carried out by multi/hyper-spectral images:

(i) spatial information contained in the position of each pixel and in the geometry of neighborhood pixels according to the location on the bitmap grid; (ii) spectral information contained in the spectral bands, in such a way that a spectrum is associated to each pixel.

Using simultaneously both kinds of information is a crucial point to get a robust and accurate method of segmentation. More precisely, one can compare the spectral information of each pixel and group them into classes. This comparison is global on the image, since each pixel is compared to all the others. However, in this point-wise comparison, the spatial relationship is usually not taken into account. Another method consists in comparing the pixels in a neighborhood leading to an approach based on a spatial comparison, but usually this spatial processing is achieved impenitently for each spectral band. With these remarks, we understand the necessity of combining both kinds of information.

State-of-the-art on segmentation of multispectral images. In the literature, there are several other methods of segmentation that can be used for segmentation of multispectral images.

Watershed based segmentation was used by [START_REF] Soille | Morphological partitioning of multispectral images[END_REF] to separate into classes the histogram of the multispectral image. Mean Shift algorithm was used to segment hyperspectral data cubes by [START_REF] Genova | Reduction and segmentation of hyperspectral data cubes[END_REF]. Moreover, various classification methods were considered to group pixels into non-connected classes. Support Vector Machines (SVM) were applied for instance to remote sensing by [START_REF] Gualtieri | Support Vector Machines for Hyperspectral Remote Sensing Classification[END_REF], [START_REF] Lennon | Classification of hyperspectral images with nonlinear filtering and support vector machines[END_REF], [START_REF] Archibald | Feature Selection and Classification of Hyperspectral Images With Support Vector Machines[END_REF]. One of the advantages of this method is not to be sensitive to dimensionality problem, also called "Hughes phenomenon" [START_REF] Hughes | On the mean accuracy of statistical pattern recognizers[END_REF][START_REF] Landgrebe | Hyperspectral Image Data Analysis[END_REF]. "K-nearest neighbours" was used as a spatio-spectral classification method by [START_REF] Marcal | Hierarchical clustering of multispectral images using combined spectral and spatial criteria[END_REF], [START_REF] Polder | Spectral imaging for measuring biochemicals in plant material[END_REF]. In [START_REF] Schmidt | Wavanglet: an efficient supervised classifier for hyperspectral images[END_REF], a classification method based on wavelets was introduced.

In order to analyse the structures in an image, [START_REF] Pesaresi | A new approach for the morphological segmentation of high resolution satellite imagery[END_REF] developed the morphological profile, based on granulometry principle [START_REF] Serra | Image Analysis and Mathematical Morphology[END_REF]. The Morphological Profile is composed of a series of openings by reconstruction of increasing sizes and of a series of closings by reconstruction (obtained by the dual operation) [START_REF] Soille | Morphological Image Analysis[END_REF]. With these operations, the size and the shape of the objects contained in the image are determined. Therefore, the MP contains only spatial information. A straightforward method to extend MP to hyperspectral image, is to build the MP for each channel. This approach is called Extended Morphological Profile (EMP). In fact, [START_REF] Benediktsson | Classification of hyperspectral data from urban areas based on extended morphological profiles[END_REF] use the first principal components obtained by a Principal Component Analysis of the hyperspectral image. Consequently, EMP contains spatial and spectral information. Finally a classifier, such as a SVM or a neural network [START_REF] Fauvel | Spectral and spatial classification of hyperspectral data using SVMs and morphological profiles[END_REF], is applied on EMP to obtain a classification in which the classes are not necessarily connected. In the present paper, we consider a segmentation as a partition of connected classes of an image following the definition given in [START_REF] Serra | A lattice approach to image segmentation[END_REF] in the framework of lattice theory.

Gradient based approaches such as the watershed segmentation cannot detect thin objects in remote sensing images. Indeed, thin objects have no interior on the image of the gradient. Consequently, the flooding procedure which defines a region starting from markers on the gradient, does not take into account the thin objects. These thin objects and other "small" structures can be obtained, on the one hand, using segmentation techniques based on connective criteria for pixel aggregation (Noyel et al. 2007b[START_REF] Soille | Constrained connectivity for hierarchical image partitioning and simplification[END_REF], or on the other hand, using the residue of an opening (known as top-hat operator) which has been extended to multivariate images in [START_REF] Angulo | Modelling and segmentation of colour images in polar representations[END_REF]. As discussed in [START_REF] Soille | Constrained connectivity for hierarchical image partitioning and simplification[END_REF], this is a limitation for some applications, e.g., extraction of road networks in aerial images.

The watershed transformation (WS) is one of the most powerful tools for segmenting images and was introduced by [START_REF] Beucher | Use of watersheds in contour detection[END_REF]. According to the flooding paradigm, the watershed lines associate a catchment basin to each minimum of the relief to flood (i.e. a greyscale image) [START_REF] Beucher | The morphological approach to segmentation: the watershed transformation[END_REF]. Typically, the relief to flood is a gradient function which defines the transitions between the regions. Using the watershed on a scalar image without any preparation leads to a strong over-segmentation (due to a large number of minima). There are two alternatives in order to get rid of the over-segmentation. The first one consists in initially determining markers for each region of interest. Then, using the homotopy modification, the only local minima of the gradient function are imposed by the region markers. The extraction of the markers, especially for generic images, is a difficult task. The second alternative involves hierarchical approaches based on non-parametric merging of catchment basins (waterfall algorithm) or based on the selection of the most significant minima according to different criteria (dynamics, area or volume extinction values) using extinction functions [START_REF] Meyer | An overview of morphological segmentation[END_REF].

We have introduced a general method to segment hyperspectral images by WS (Noyel et al. 2007a). Several multivariate gradients were studied. The markers are obtained from a previous spectral classification.

Although this approach is powerful, it is a deterministic process which tends to build irregulars contours and can produce an over-segmentation. The stochastic WS was proposed in order to regularize and to produce more significant contours [START_REF] Angulo | Stochastic watershed segmentation[END_REF]. The initial framework was then extended to multispectral images (Noyel et al. 2007c, 2008a[START_REF] Noyel | Filtrage, réduction de dimension, classification et segmentation morphologique hyperspectrale[END_REF].

Aim of the paper. Based on the experience of our previous works on multivariate image segmentation, the aim of this study is to present a complete chain for an automatic spatio-spectral segmentation of multispectral images, and to illustrate the method with some examples of remote sensing images.

Let us start by presenting a general paradigm of WS-based segmentation of multispectral images (fig. 1). This segmentation requires two basic ingredients:

(i) some markers for the regions of interest mrk;

(ii) a relief to flood which describes the "energy" of the frontiers between the regions g.

The markers can be chosen interactively by a user or automatically by means of a morphological criteria ξ N , or with the classes of a previous spectral classification.

The relief to flood is a scalar function (i.e., a greyscale image). For the deterministic WS, it is usually a gradient (in fact its norm), or a distance function. For the stochastic WS, the function to flood is a probability density function (pdf) of the contours appearing in the image. After extracting the markers, they are imposed as sources of the relief to flood and the WS is computed. The results are noted W S(g, mrk) or W S(g, ξ N ).

The segmentation framework for remote sensing images introduced in this paper is coherent with the general paradigm of WS-based segmentation. More precisely, the method is based on classification-driven stochastic watershed by Monte Carlo simulations, and as we will show, it gives more regular and reliable contours than standard watershed.

The approach is decomposed into several sequential steps. Each step was partially considered in our previous works, where our methods were compared with more standard ones. However, the main objective of this paper is to introduce the global pipeline, presenting exclusively for each step the algorithm leading to the best performance.

Firstly in section 3., a dimensionality reduction stage is performed using Factor Correspondence Analysis (FCA) [START_REF] Benzécri | L'Analyse Des Données[END_REF]. In this context, a new way to select the factor axes (i.e. the eigenvectors) according to their spatial information is introduced. Then a spectral classification on the factor space, described in section 5., produces a spectral pre-segmentation of the image. Subsequently, a probability density function (pdf) of contours containing spatial and spectral information is estimated by simulation using a stochastic WS approach driven by the spectral classification. In section 4., the standard stochastic WS is reminded, and then, in section 6., the algorithm for the construction of the classification-driven stochastic WS is discussed. It also includes how the pdf of contours is finally segmented by a watershed controlled by markers coming from a regularization of the initial classification.

The paper is completed with section 2., which fixes the notations used in this paper and which presents the image data considered in the experiments. Finally, in section 7., the conclusions and perspectives of the paper are discussed.

Notations

In a formal way, each pixel of a multispectral image is a vector with values in wavelength. To each wavelength corresponds an image in two dimensions called channel. The number of channels depends on the nature of the specific problem under studies (satellite imaging, spectroscopic images, temporal series, etc.). Let

f λ : E → T L x → f λ (x) = (f λ1 (x), f λ2 (x), . . . , f λL (x)) (1)
denote an hyperspectral image, where:

• E ⊂ R 2 , T ⊂ R and T L = T × T × . . . × T • x = x i / i ∈ {1
, 2, . . . , P } is the spatial coordinate of a vector pixel f λ (x i ) (P is the number of pixels of E) • f λj / j ∈ {1, 2, . . . , L} is a channel (L is the number of channels) • f λj (x i ) is the value of vector pixel f λ (x i ) on channel f λj . Figure 2 gives an example of a five band satellite simulated image PLEIADES, acquired by the CNES (Centre National d'Etudes Spatiales, the French space agency) and provided by [START_REF] Flouzat | Spatial and spectral segmentation of satellite remote sensing imagery using processing graphs by mathematical morphology[END_REF]. Its channels are the following: f λ1 blue, f λ2 green, f λ3 red, f λ4 near infrared and f λ5 panchromatic. The panchromatic channel, initially 1460×1460 pixels with a resolution of 0.70 meters, was resized to 365×365 pixels. Therefore, the resolution is 2.80 meters in an image of 365×365×5 pixels. In order to represent a multispectral image in a synthetic way, we created a synthetic RGB image using channels f λ3 red, f λ2 green and f λ1 blue.

The downsampling method consists in, by averaging, keeping only one pixel on four pixels in both directions of the plane. The coordinates (i, j) of the pixels are ordered as the standard ordering used for matrices. Moreover, we also notice that even if the downsampled panchromatic image has values strongly correlated to average intensity of the red, blue and green channels, the spatial resolution of the panchromatic is better than the chromatic bands, and consequently there is a better resolution of the contours which is important for the aim of image segmentation. Indeed, to deal with the redundancy, the FCA (Factor Correspondence Analysis) transforms the original channels into uncorrelated factor channels (according to the chi-squared metric).

(a) f λ1 (b) f λ2 (c) f λ3 (d) f λ4 (e) f λ5 (f ) RGB

Dimensionality reduction

As the number of channels, in a multispectral image, is usually important, the spectrum information contained in these channels is redundant. Therefore, in order to avoid the Hughes phenomenon [START_REF] Hughes | On the mean accuracy of statistical pattern recognizers[END_REF][START_REF] Landgrebe | Hyperspectral Image Data Analysis[END_REF]) associated to the spectral dimensionality problem as well as to reduce the computational time, it is necessary to reduce the number of channels. Consequently, a data reduction is performed using Factor Correspondence Analysis (FCA) [START_REF] Benzécri | L'Analyse Des Données[END_REF]. We prefer using a FCA instead of a Principal Component Analysis (PCA) or a 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], because image values are positive and the spectral channels can be considered as probability distributions. [START_REF] Green | A transformation for ordering multispectral data in terms of image quality with implications for noise removal[END_REF] have demonstrated that PCA and MNF are equivalent in the case of an uncorrelated noise with equal variance in all channels. This explains why PCA is often used to remove noise on hyperspectral images in which the noise is almost uncorrelated with the same variance on all channels.

As for PCA, from selected factorial FCA axes the image can be partially reconstructed. The metric used in FCA is the chi-squared, which is adapted to probability laws and normalized by channels weights. FCA can be seen as a transformation going from image space to a factorial space. In the factorial space, the coordinates of the pixel vector, on each factorial axis, are called pixel factors. The pixel factors can be considered as another multispectral image whose channels correspond to the factorial axes:

ζ : T L → T K / K < L f λ (x) → c f α (x) = c f α1 (x), . . . , c f αK (x) (2)
A limited number K, with K < L, of factorial axes is usually chosen. Therefore FCA can be seen as a projection of the initial vector pixels in a factor space with a lower dimension. Moreover, as for PCA, FCA reduces the spectral noise on multispectral images [START_REF] Green | A transformation for ordering multispectral data in terms of image quality with implications for noise removal[END_REF], Noyel et al. 2007a[START_REF] Noyel | Filtrage, réduction de dimension, classification et segmentation morphologique hyperspectrale[END_REF]. Figure 3 shows FCA factor pixels c f α of image "Roujan". Consequently, we have two spaces for multivariate segmentation: the multispectral image space (MIS) and the factor image space (FIS).

A common problem in remote sensing application is the registering of the various spectral images. If the registering step is not perfect, the spatial shifts introduce in the dimensionality reduction a FCA factor containing this information; see for instance, in figure 3, the FCA factor image c f α2 . This image looks like a laplacian image, but for the purpose of segmentation, it should be considered as a factor image containing spatial noise. We only show below how to detect and remove these "noisy" factor images.

Another approach to reduce the number of channels consists in modelling the spectrum of each vector-pixel f λ (x i ) and creating a new multispectral image composed of the parameters of the model (Noyel et al. 2007a); however, for generic multispectral remote sensing images, this kind of modelling is quite difficult. As shown by [START_REF] Benzécri | L'Analyse Des Données[END_REF] and by [START_REF] Green | A transformation for ordering multispectral data in terms of image quality with implications for noise removal[END_REF], some factor pixels on factor axes contain mainly noise, and others signal information. In order to select relevant axes containing information (i.e. signal against noise), we introduce a new method based on a measurement of the signal to noise ratio SN R of every factor image. For each factorial axis α k , the centered spatial covariance is computed by a 2D FFT (Fast Fourier Transform) on the pixel factors:

g α1 (h) γg α1 (h) g α2 (h) γg α2 (h)
g αk (h) = E[c f αk (x)c f αk (x + h)], with c f αk (x) = c f αk (x) -E[c f αk (x)]
, where E is the mathematical expectation (here the statistical expectation).

The covariance peak, at the origin, contains the sum of the signal variance and the noise variance of the image. Then, the signal variance is estimated by the maximum (i.e. value at the origin) of the covariance g. Based on the property of the morphological opening which removes intensity peaks smaller than the used structuring element, we apply a morphological opening γ on the covariance image with a unitary centered structuring element (square of 3 × 3 pixels) in order to remove the peak of signal associated to the noise variance: V ar(signal) = γg αk (0). The noise variance is therefore given by the residue of the opening of the covariance at the origin: V ar(noise) = g αk (0) -γg αk (0) (fig. 4). Hence, the signal to noise ratio is defined for a factor axis α k as:

SN R αk = V ar(signal) V ar(noise) = γg αk (0) g αk (0) -γg αk (0)
(3)

In the current example, by analysis of the factor pixels and their signal to noise ratios, we observe that they are higher for axes 1, c f α1 , and 3, c f α3 , than for axes 2, c f α2 , and 4, c f α4 . Therefore, the axes 1 and 3 are retained. For an automatic selection of factor axes, the axes with a SNR lower than 1 are rejected.

In figure 5, we notice that axis 3, which is selected, has a lower inertia than axis 2, which is rejected. Therefore, SNR analysis makes it possible to describe relevant signal in dimensionality reduction, more than inertia-like criterion.

4. Introduction to stochastic WS [START_REF] Angulo | Stochastic watershed segmentation[END_REF] defined a new method of stochastic WS for greyscale and color images. This method was extended to hyperspectral images by Noyel et al. (2007c). In addition, the segmentation of multispectral images by classical (deterministic) WS was presented in (Noyel et al. 2007a). Several improvements were made by Noyel et al. (2008a), [START_REF] Noyel | Filtrage, réduction de dimension, classification et segmentation morphologique hyperspectrale[END_REF].

Principle of the stochastic WS

One of the main artefacts of the classical watershed is that small regions strongly depend on the position of the markers, or on the volume (i.e. the integral of the grey levels) of the catchment basin, associated to their minima. In fact, there are two kinds of contours associated to the watershed of a gradient: first order contours, which correspond to significant regions and which are relatively independent from markers; second order contours, associated to "small", "low" contrasted and textured regions, which depend strongly on the location of markers. Stochastic watershed aims at enhancing the first order contours from a sampling effect, to improve the result of the watershed segmentation.

Let us consider {mrk i (x)} M i=1 a series of M realisations of N uniform or regionalized random germs. Each of these binary images is considered as the marker for a watershed segmentation of a scalar gradient or a vector gradient. Therefore, a series of segmentations is obtained, i.e. {sg mrk i (x)} M i=1 . Starting from the M realisations of contours, the probability density function pdf (x) of contours is computed by the Parzen window method. The kernel density estimation by Parzen window [START_REF] Duda | Pattern Classification and Scene Analysis[END_REF]) is a way to estimate the probability density function (pdf) of a random variable. Let x 1 , x 2 , • • • , x M ∈ R n be M samples of a random variable, the kernel density approximation of its pdf is:

f h (x) = 1 N h N i=1 K x-xi h ,
where K(x) is some kernel and the bandwidth h a smoothing parameter. Usually, K(x) is taken to be a Gaussian function, G σ , with mean zero and variance σ 2 .

The smoothing effect of the Gaussian convolution kernel (typically σ = 3 working on contours of one pixel width) is important to obtain a function where near contours, such as textured regions or associated to small regions, are added together. In other words, the WS lines with a very low probability, which correspond to non significant boundaries, are filtered out.

The pdf of contours could be thresholded to obtain the most prominent contours. However, to obtain closed contours, the pdf image can be also segmented, as we consider here, using a watershed segmentation.

Influence of parameters

The stochastic WS needs two parameters:

(i) M realisations of germs. The method is almost independent on M if it is large enough. Practically, the convergence is ensured for M in the range 20 -50. We propose to use M equal to 100. (ii) N germs (or markers): In WS segmentation the number of regions obtained is equal to the number of markers. In the case of the stochastic WS, if N is small, a segmentation in large regions is privileged; if N is too large, the over-segmentation of sg mrk i leads to a very smooth pdf , which looses its properties to select the R regions. Therefore, the stochastic WS mainly depends on N , but is linked to the number of regions R which should be finally selected. As it was shown in [START_REF] Angulo | Stochastic watershed segmentation[END_REF], it is straightforward to use N > R.

As it is described below, in section 6., in this study the number of markers of each realization used in the conditioned stochastic WS will be equal to the number of connected classes of the classification κ different of the void class.

Generation of markers from a classification

To segment the pdf of contours by WS, it is necessary to have markers for the spectral objects of interests. The latter are obtained by a spectral classification. This classification step groups the pixels into classes with similar spectra. In fact, each pixel is compared to all the others, leading to a global comparison of the image spectra. We have tested two kinds of methods for spectral classification on multispectral images: unsupervised methods and supervised methods.

In the case of remote sensing, we show results obtained by an unsupervised method such as "clara" [START_REF] Kaufman | Finding Groups in Data[END_REF]. Unsupervised and supervised methods in remote sensing using mathematical morphology as well as the "clara" method are also investigated in [START_REF] Epifanio | Morphological texture features for unsupervised and supervised segmentations of natural landscapes[END_REF]. For other contexts different from remote sensing, we tested other unsupervised methods such as k-means and supervised methods such as Linear Discriminant Analysis (Noyel et al. 2008b[START_REF] Noyel | Filtrage, réduction de dimension, classification et segmentation morphologique hyperspectrale[END_REF]. One of the strong points of supervised approaches is to integrate prior information about spectrum classes into the classification. However, in generic remote sensing segmentation problems, we are not looking for a special kind of spectrum. Therefore, we prefer, in the presentation of the current paper, to use an unsupervised classification with a unique parameter, namely the number of classes.

Hence, after the data reduction and the spectral filtering stage by FCA, a spectral classification by "clara" is performed on the factor space. We have to stress that a correct classification, used later for constructing the pdf, requires factor pixels without noise. The classification algorithm uses an Euclidean distance which is coherent with the metric of the factor space (Noyel et al. 2007a). The unique parameter is the number of classes, which is chosen in order to get a good separation of the classes. Here, for this example, 3 classes are retained (fig. 6).

Multispectral segmentation by stochastic watershed

The previously obtained spectral classification of pixels leads to a spatial partition of the image on regions which are spectrally homogeneous. This initial partition is denoted κ. However, they do not define a partition into spatial connected classes: various connected components belong to the same spectral class. In order to address this issue, a segmentation stage is needed which introduces spatial information (i.e., regional information). In this section, after presenting an algorithm to pre-process the previous classification κ, the extension of stochastic WS to multispectral images is explained; then we show a way to constrain the pdf by the previous spectral classification.

Pre-processing of markers coming from the classification

A pre-processing of the connected classes of the classification is necessary for two reasons:

(i) to give necessary degrees of freedom to the final WS of the pdf. Actually, if we use the connected classes of the classification κ as markers of the flooding process, then the WS is completely defined by the limit of the classes after labelling them into connected classes. Therefore we propose to process the initial partition using an anti-extensive transformation such as an erosion [START_REF] Serra | Image Analysis and Mathematical Morphology[END_REF][START_REF] Soille | Morphological Image Analysis[END_REF]. For instance, using a structuring element (SE) of size 5 × 5 pixels each spatial class is reduced and the smallest classes disappear. We consider that classes corresponding to "noise" are totally removed. The smallest classes could also be removed by means of an area opening [START_REF] Soille | Morphological Image Analysis[END_REF]. Therefore, after processing the initial partition, a partial partition is obtained. There are several alternatives to deal with the problem that the partial partition is not a partition of the support space; for instance adding to the partial partition the singletons outside its support. We prefer however the following solution. We introduce a particular label, the "void class", in such a way that after processing the initial partition, covering the whole space, a new partition is also obtained. In the new partition, the existing classes are modified (reduced) and in addition, a new class is introduced, the void class. If we want to keep all the classes, for instance if the aim is to segment thin objects as a network of roads, we could use a homotopic thinning instead of an erosion. (ii) to fill the holes inside the largest classes with an extensive transformation such as a closing by reconstruction with a SE of size for instance 3 × 3 pixels [START_REF] Serra | Image Analysis and Mathematical Morphology[END_REF][START_REF] Soille | Morphological Image Analysis[END_REF].

The transformations are applied to each class of the initial spectral classification. More precisely to each spectral class, κ(n) with n = 1, . . . , N , of the previous classification κ = ∪ n κ(n), an index function (i.e. a binary image) is associated:

h κ(n) (x) = 1, if x ∈ κ(n); 0, otherwise. (4)
Then, the sequence of the anti-extensive and extensive transformations are applied to each index function, h κ(n) for all n = 1, . . . , N , in order to obtain the transformed index function h κ(n) . Then the supremum (union) of the index functions is computed, followed by a standard labelling to obtain the connected components which correspond to markers for the segmentation, except for the void class which is not a marker. The complete transform is noted Υ, and it results in a new partition κ = Υ(κ) (fig. 6). The partition κ will then be strongly used for the spatial segmentation by stochastic WS.

Extension of the stochastic WS to multivariate or hyperspectral images

The extension of the stochastic WS was first presented in Noyel et al. (2007c). The key points are recalled in the current paper.

6.2.1

Spectral distances and gradient. In order to segment images according to watershed-based paradigms, a gradient is needed. A gradient image, in fact its norm, is a scalar function with values in the reduced interval [0, 1] (after normalisation), i.e. (x) : E → [0, 1]. In order to define a gradient, two approaches are considered: the standard symmetric morphological gradient on each marginal channel and a metric-based vectorial gradient on all channels.

The morphological gradient is defined for scalar images f as the difference be-tween dilation and erosion by a unit structuring element B, i.e.,

(f λj (x)) = δ B (f λj (x)) -ε B (f λj (x)) = ∨[f λj (y), y ∈ B(x)] -∧[f λj (y), y ∈ B(x)] .
As it was shown in [START_REF] Hanbury | Morphological operators on the unit circle[END_REF], the morphological gradient can be generalised to vectorial image by using the following relation:

{f (x) -∨[f (y), y ∈ B(x), y = x]} = ∧{f (x) -f (y), y ∈ B(x), y = x}
and the relation obtained by the inversion of the suprema and infima. The morphological gradient can be written in a form that is only composed of increments computed in the neighborhood B centered at point x:

(f )(x) = ∨[f (x) -f (y), y ∈ B(x), y = x] -∧[f (x) -f (y), y ∈ B(x), y = x] .
To extend this relation to multivariate functions, it is sufficient to replace the increment [f (x)-f (y)] by a distance between vector pixels d(f λ (x), f λ (y)) to obtain the following metric-based gradient:

d f λ (x) = ∨[d(f λ (x), f λ (y)), y ∈ B(x), y = x] -∧[d(f λ (x), f λ (y)), y ∈ B(x), y = x] .
Various metric distances, useful for multispectral images, are available for this gradient such as:

• the Euclidean distance:

d E (f λ (x), f λ (y)) = L j=1 (f λj (x) -f λj (y)) 2 ,

and

• the Chi-squared distance:

d χ 2 (f λ (x i ), f λ (x i )) = L j=1 S f .λj f λj (x i ) f xi. - f λj (x i ) f x i . 2 with f .λj = P i=1 f λj (x i ), f xi. = L j=1 f λj (x i ) and S = L j=1 P i=1 f λj (x i
). An important point is to choose an appropriate distance depending on the space used for image representation: Chi-squared distance is adapted to MIS and Euclidean distance to FIS. More details on multivariate gradients are given by Noyel et al. (2007a). Another example of a multivariate gradient is given by Scheunders (2002).

6.2.2

Probability density function for multispectral images. We studied two ways to extend the pdf gradient to multispectral images: (i) the first one is a marginal approach (i.e. channel by channel) called marginal pdf mpdf (alg. 1 and fig. 7) (ii) the second one is a vectorial approach (i.e. vector pixel by vector pixel) called vectorial pdf vpdf (alg. 2 and fig. 8). ()

) ( , ), ( ), ( ) ( 2 1 x f x f x f x L λ λ λ λ K = f R regions … ) ( 1 1 x mrk ) ( 1 x mrk M … ) ( 2 1 x mrk ) ( 2 x mrk M … ) ( 1 x mrk L ) (x mrk L M … ) ( 1 1 x sg ) ( 1 x sg M ∑ = * = M i i G x sg M x pdf 1 1 1 ) ( 1 ) ( σ PDF (Parzen method) … ) ( 2 1 x sg ) ( 2 x sg M PDF (Parzen method) … ) ( 1 x sg L ) (x sg L M ∑ = * = M i L i L G x sg M x pdf 1 ) ( 1 ) ( σ space image in / 1 , L w j j = ∀
vpdf (x) = 1 M × L M ×L i=1 sg i (x) * G σ .
(6)

PDF (Parzen method)

Vectorial gradient

Random germs

Marker-controlled WS

mrk 1 (x) mrk 2 (x) mrk M×L (x) sg 1 mrk (x) sg 2 mrk (x) sg MxL mrk (x) sg R-vol (vpdf)(x) ∑ × = * × = L M i i G x sg L M x vpdf 1 ) ( 1 ) ( σ (
) The probabilistic gradient was also defined in [START_REF] Angulo | Stochastic watershed segmentation[END_REF] to ponder the enhancement of the largest regions by the introduction of smallest regions. It is defined as prob = mpdf + d : after normalization in [0, 1] of the weighted marginal pdf mpdf and the metric-based gradient d .

) ( , ), ( ), ( ) ( 2 1 x f x f x f x L λ λ λ λ K = f Multispectral image (vectorial)

R regions

In order to obtain a partition from the mpdf , vpdf or the gradient prob , these probabilistic functions can be segmented for instance by a hierarchical WS with a volume criterion (see the examples of fig. 9), as studied in Noyel et al. (2007c). In such a case, the goal is not to find all the regions. In fact, the stochastic WS addresses the problem of segmentation of an image in few pertinent regions, according to a combined criterion of contrast and size. In the present study, as we discuss below, the segmentation of the pdf is obtained from markers of the processed classification. In Noyel et al. (2007c) was also shown that marginal pdf and vectorial pdf give similar results in multispectral image space (MIS) and in factor image space (FIS), as we can also observe in comparison of fig. 9. Consequently, in what follows, we choose to work on the MIS to show that our method is not limited by the number of channels.

(a) χ 2 (f λ ) (b) mpdf (f λ ) (c) vpdf (f λ ) (d) prob (f λ ) (e) seg vol ( χ 2 (f λ ), R) (f ) seg vol (mpdf (f λ ), R) (g) seg vol (vpdf (f λ ), R) (h) seg vol ( prob (f λ ), R)

Conditioning the germs of the pdf by a previous classification

The pdf of contours with uniform random germs is constructed without any prior information about the spatial/spectral distribution of the image. In this part, we introduce spectral information by conditioning the germs by the previous transformed classification κ. Therefore this obtained pdf contains spatio-spectral information. It is possible to use point germs or random ball germs whose location is conditioned by the classification. An exhaustive study is presented by [START_REF] Noyel | Filtrage, réduction de dimension, classification et segmentation morphologique hyperspectrale[END_REF], Noyel et al. (2008c). In the sequel we present random ball germs regionalized by a classification where each connected class may be hit one time, mrk κ-b i (x). Let us explain the procedure. The transformed classification κ is composed of connected classes, κ = ∪ k C k with C k ∩C k = ∅, for k = k . The void class is written C 0 . Then the random germs are drawn conditionally to the connected components C k of the filtered classification κ. To do this, the following rejection method is used: the random point germs are uniformly distributed. If a point germ m falls inside a connected component C k of minimal area S and not yet marked, then it is kept, otherwise it is rejected. Therefore not all the germs are kept. These point germs are called random point germs regionalized by the classification κ. However, the regionalized point germs are sampling all the classes independently of their prior estimate of class size/shape, given by the classification. In order to tackle with this limitation, we propose to use random balls as germs.

The centers of the balls are the random point germs and the radii r are uniformly distributed between 0 and a maximum radius Rmax: U[1, Rmax]. Only the intersection between the ball B(m, r) and the connected component C p is kept as a germ. These balls are called random balls germs regionalized by the classification κ and noted mrk κ-b i (x). The algorithm 3 sketches the process. We notice that N is the number of random germs to be generated. The effective number of implanted germs is less than N .

We use the marginal pdf to show our results. Some random ball germs regionalized by a classification, {mrk κ-b,j i (x)} j=1 i=1...5 , their associated realisations of contours, {sg j i (x)} j=1 i=1...5 , and the marginal pdf computed in MIS, mpdf κ-b (f λ ), are presented in figure 10. Set the class C k to marked 8: A comparison between the marginal pdf with uniform random point germs, mpdf pt (f λ ) (noted before mpdf (f λ )) and the marginal pdf with random ball germs mpdf κ-b (f λ ) is shown in figure 11. These pdf are segmented by hierarchical WS with a volume criterion, sg R-vol (mpdf pt (f λ )) and sg R-vol (mpdf κ-b (f λ )). In this particular example, the advantage of the pdf obtained with the conditioned random balls is obvious with respect to the pdf from uniform random point germs. Moreover, in order to compare the results of segmentation of example of figure 11, we must notice that in both images only the 8 most important regions have been segmented. If the purpose is to segment all the objects, including the smallest ones, we only need to increase the number of desired volumic regions.

end if 9: end for i = 1 i = 2 i = 3 i = 4 i = 5 (a) mrk κ-b,1 1 (x) (b) mrk κ-b,1 2 (x) (c) mrk κ-b,1 3 (x) (d) mrk κ-b,1 4 (x) (e) mrk κ-b,1 5 (x) (f ) sg 1 1 (x) (g) sg 1 2 (x) (h) sg 1 3 (x) (i) sg 1 4 (x) (j) sg 1 5 (x) (k) mpdf κ-b (f λ )
(a) mpdf pt (f λ ) (b) sg R-vol (mpdf pt (f λ )) (c) mpdf κ-b (f λ ) (d) sg R-vol (mpdf κ-b (f λ ))

Segmentation of the pdf

After having built the pdf, a segmentation stage is required to obtain the contours of the image. Two possibilities were tested:

(i) a marker-controlled WS using as markers the transformed classification sg mrk (mpdf κ-b (f λ ), κ) (Noyel et al. 2008a). In this case, no more parameter is required. (ii) a hierarchical WS with a volume criterion sg R-vol (mpdf κ-b (f λ )) (Noyel et al. 2007c). Therefore, it is necessary to define an additional parameter: the number of regions with the largest volume.

In order to show the robustness of our approach we used the same method to build the pdf conditioned by the classification on similar images "Roujan", "Roujan 0 2" and "Roujan 0 9". As for image "Roujan", factor axes c f α1 and c f α3 are kept because th Signal to Noise Ratio (SNR) is less than thresholded value (1). Inertias and SNR for factor axes of images "Roujan X X" are presented in the tables 1 and 2.

A classification "clara" is made on the factor pixels of these two axes. The marginal pdf, mpdf κ-b (f λ ), is produced as explained for image "Roujan".

Table 1. Part of inertia for factor axes of images "Roujan X X"

Image c f α 1 c f α 2 c f α 3 c f α 4
"Roujan" 84.1 % 8.7 % 6.2 % 1 % "Roujan 0 2" 75.6 % 13.7 % 9.2 % 1.5 % "Roujan 1 9" 77.5 % 12.1 % 9 % 1.4 % Table 2. SNR for factor axes of images "Roujan X X"

Image c f α 1 c f α 2 c f α 3 c f α 4
"Roujan" 5.10 0.03 2.84 0.47 "Roujan 0 2" 2.66 0.01 2.53 0.32 "Roujan 1 9" 3.15 0 2.90 0.51

The parameters for the stochastic WS controlled by markers sg mrk (mpdf κ-b (f λ ), κ) are:

• the number of classes for the classification κ: Q = 3. It is the only important parameter; • the maximum number of random ball germs: N = 50, which must be high enough and which must be of the the same order as the number of regions in the image; • the number of realisations for each channel: M = 100, which always has the same value; • the minimum area S = 10 pixels for the connected classes of the classification, generally the same; • the parameters for pre-processing the initial spectral classification are an erosion using a 5 × 5 square and a closing by reconstruction using as marker a dilation using a 3 × 3 square, generally the same; • the maximum radius of the random ball Rmax = 30, generally the same.

The final segmentations with marker-controlled WS, sg mrk (mpdf κ-b (f λ ), κ), and hierarchical WS with a volume criterion, sg R-vol (mpdf κ-b (f λ )), are shown on figure 12. We notice that the number of regions, resulting from the segmentation, strongly depends on the image, whereas the number of spectral classes is the same for similar images. In fact, the number of regions for volumic WS, sg vol , must be chosen according to the spatial diversity of the considered image, and cannot be easily fixed "a priori" (50 here). The number of regions depends on the size and the complexity of the image, while the number of spectral classes depends on the spectral content, which for a specific domain can be relatively constant. Consequently, it is more relevant to segment the pdf with markers coming also from the classification sg mrk . Indeed, only one parameter is needed for marker-controlled WS sg mrk (mpdf κ-b (f λ ), κ): the number of classes Q in the classification κ. Moreover, in that sense, the approach based on the number of classes Q produces a more robust segmentation than the approach based on the number of regions R.

In order to show that the stochastic WS, conditioned by a previous classification, produces contours which are more regular and robust than deterministic (standard) WS, we compare the results of a deterministic approach to a stochastic one:

(i) For both cases, an unsupervised classification κ ("clara") is processed in factor space FIS and transformed κ = Υ(κ). (ii)

• For the deterministic approach, a chi-squared metric based gradient χ 2 (f λ ) is computed in image space MIS, as a function to flood. • For the stochastic approach, a marginal probability density function mpdf κ-b (f λ ), with regionalized random ball germs conditioned by the classification κ, is processed, in image space MIS, as a function to flood. (iii) In both cases, the flooding function is segmented by a watershed (WS) using as sources of flooding the markers from the transformed classification κ.

The parameters for the stochastic WS controlled by markers are the same as in part 6.4. In figure 13, are given the results of segmentation by deterministic and stochastic WS for image "Roujan" with markers coming from the transformed classification κ. We observe that the contours are smoother and follow more precisely the main limits of the regions for the stochastic approach than for the deterministic one. The same observation is made on the image "Port de Bouc" classified in Q = 7 classes. (fig. 14). 6.5 General method of multispectral segmentation by stochastic WS After all these observations and analysis, we can summarize now the algorithmic pipeline of the general method of segmentation of multispectral image by stochastic WS. The framework is presented on figure 15. The different steps are as follows:

c) χ 2 (f λ ) (d) sg mrk ( χ 2 (f λ ), κ) (e) zoom sg mrk ( χ 2 (f λ ), κ) (f ) mpdf κ-b (f λ ) (g) sg mrk (mpdf κ-b (f λ ), κ) (h) zoom sg mrk (mpdf κ-b (f λ ), κ) (
c) χ 2 (f λ ) (d) sg mrk ( χ 2 (f λ ), κ) (e) zoom sg mrk ( χ 2 (f λ ), κ) (f ) mpdf κ-b (f λ ) (g) sg mrk (mpdf κ-b (f λ ), κ) (h) zoom sg mrk (mpdf κ-b (f λ ), κ) (
• a FCA to reduce the spectral dimension of the image f λ and eventually to filter the noise. This transformation creates a factor image c f α

• the previous step requires the selection of factorial axes by the signal to noise ratio (SNR) of the factor pixels • a spectral classification to group pixels into homogeneous classes (global pointwise comparison) • a morphological transform Υ of each classes. The image transformed is noted κ = Υ(κ) • the building of a pdf of contours that can be marginal mpdf or vectorial vpdf starting from a gradient on the image space f λ (MIS) or factor space c f α (FIS) with random-balls germs regionalized by the spectral classification • a segmentation of the pdf by a WS controlled by markers coming from the classes of the transformed classification κ Some of the steps presented in the flowchart of figure 15 have not been discussed in this paper. These steps are very useful for other applications of multivariate image segmentation; the interested reader can find more details in Noyel et al. (2008b). The part outlined with a dotted line is not presented in this paper.

f   x= f  1 , , f  L  FCA ζ c  f x =c  1 f , , c  K f    

Conclusion

A novel method to segment multispectral images is presented in this paper. First, a dimensionality reduction by Factor Correspondence Analysis is performed. It reduces the spectral noise and preserves the spatial information. A new way to select factor axes is presented. It uses the spatial information carried by the factor pixels (i.e. signal to noise ratio) and not on only the statistical information (i.e. inertia). Then, a classification stage ensures a global pixel-wise comparison and groups pixels into classes with similar spectra. A morphological transformation is subsequently applied to obtain markers for the final segmentation. Then, a probability density function of contours conditioned by the previous classification is built with random ball-germs. This pdf contains spatial and spectral information and is estimated by Monte-Carlo simulations. During the pdf estimation process regional information between pixels is compared. Finally, the pdf is segmented by a WS controlled by markers coming from the transformed classification. Moreover, this spatio-spectral method of segmentation needs few and well controlled parameters:

• the threshold of the signal to noise ratio for the Factor Correspondence Analysis (dimensionality reduction stage); this can be selected according to the knowledge of the amount of noise in the images; • the number of classes of the classification which corresponds to the number of different kinds of spectra in the image. This parameter is more robust than the number of regions with the largest volume in the standard hierarchical watershed segmentation; • eventually, a size criterion to remove small regions which must be excluded of the segmentation; • the number of random germs (linked to the number of regions) N ;

• the number of realisations (fixed) M = 100;

• the maximum radius of the ball (related to the surface area of interest zones)

Rmax.

It is shown, that in the case of segmentation of similar images these parameters can be easily fixed.

The approach is perfectly valid for VHR images: the only requirement is to start from an initial classification representing all the most significant regions, and of course, the regions must have a minimal thickness to define the notion of inner maker for the watershed segmentation.

This general method of segmentation is not limited to the field of remote sensing but can also be applied on medical images (Noyel et al. 2008b), microscopy images, thermal images (Noyel et al. 2007a), time series, multivariate series, etc.

We are now thinking about two interesting perspectives: to introduce more prior information on the construction of the pdf of contours and to use a more advanced multivariate kernel during the estimation process of the pdf.

Figure 1 .

 1 Figure 1. General framework of multispectral image segmentation

Figure 2 .

 2 Figure 2. Channels of multispectral image f λ "Roujan" ( c CNES): (a) f λ 1 blue, (b) f λ 2 green, (c) f λ 3 red, (d) f λ 4 near infra-red, (e) f λ 5 panchromatic, (f ) synthetic RGB representation.

Figure 3 .

 3 Figure 3. FCA factors pixels c f α of image "Roujan" on axes (a) 1, (b) 2, (c) 3 and (d) 4 with respective inertias 84.1 %, 8.7 %, 6.2 %, 1 %.

Figure 4 .

 4 Figure 4. 2D covariances for factor pixels c f α 1 and c f α 2 of image "Roujan" before (left) and after opening (right).

Figure 5 .

 5 Figure 5. Inertia (%) and SNR for factor pixels of image "Roujan".

Figure 6 .

 6 Figure 6. First line: (a) Synthetic RGB image "Roujan", (b) classification κ in 3 classes in factor space formed by c f α 1 and c f α 3 and (c) transformed classification κ = Υ(κ). Second line: (d) Synthetic RGB image "Port de Bouc", (e) classification κ in 7 classes in factor space formed by c f α 1 and c f α 2 and (f ) transformed classification κ = Υ(κ). The void class is in black on the images κ. The way to compute the transformed classification is explained in section 6.1. The colours correspond to the labels of the transformed classification obtained after a new labellisation.

Figure 7 .

 7 Figure 7. General framework stochastic WS on marginal pdf mpdf for multispectral images

Figure 8 .

 8 Figure 8. General framework stochastic WS on vectorial pdf vpdf for multispectral images

Figure 9 .

 9 Figure 9. Top: (a) Gradients and (b, c, d) pdf; bottom: (e, f, g, h) associated WS segmentations, with a volume criterion, on image "Roujan", in MIS, with N = 50 points , M = 100 realisations, R = 50 regions with the largest volume.

Algorithm 3

 3 Random ball germs regionalized by a classification (each connected class may be hit one time) mrk κ-b i (x) 1: Given N the number of drawn germs m 2: Set the background class and the void class C 0 to marked 3: for all drawn germs m from 1 to N do 4: if C k , such as m ∈ C k , is not marked then 5: r = U[1, Rmax] 6: Keep as a germ B(m, r) ∩ C k 7:

Figure 10 .

 10 Figure 10. (a, b, c, d, e) A few realisations of regionalized random ball germs {mrk κ-b,j i (x)} j=1 i=1...5 (top) by the classification κ and (f, g, h, i, j) associated contours of watershed {sg j i (x)} j=1 i=1...5 (bottom). i is the index of realisations and j the number of the considered channel. The parameters are the maximum radius Rmax = 30 pixels and the maximum number of germs N = 50. (k) The marginal pdf with regionalized random ball germs mpdf κ-b (f λ ) is computed with M = 100 realisations.

Figure 11 .

 11 Figure 11. Comparison between (a) the pdf mpdf pt (f λ ) and (b) the associated stochastic WS with uniform random point germs with a volume criterion sg R-vol (mpdf pt (f λ )) to (c) the pdf mpdf κ-b (f λ ) and (d) the associated stochastic WS with regionalized random ball germs with a volume criterion sg R-vol (mpdf κ-b (f λ )). Only R = 8 regions with the largest volume are retained. For visualization, the contours of the segmentations are dilated by a SE of size 3 × 3 pixels.

  Figure 12. Comparison between segmentations sg mrk (mpdf κ-b (f λ ), κ) (b, e, h) with markers coming from the filtered classification κ (a, d, g) and segmentations sg R-vol (mpdf κ-b (f λ )) (c, f, i) by stochastic WS with a prior given number of regions R = 50. The results are given for images: "Roujan", "Roujan 0 2", "Roujan 1 9".

Figure 13 .

 13 Figure 13. (First line) (a) Classification in 3 classes κ and (b) markers coming from the classification κ on image "Roujan". (Second line) (d, e) deterministic approach sg mrk ( χ 2 (f λ ), κ) on (c) the chi-squared metric based gradient χ 2 (f λ ). (Third line) (g, h) Stochastic approach sg mrk (mpdf κ-b (f λ ), κ) on (f ) the marginal pdf mpdf κ-b (f λ ). The parameters are Q = 3 classes, N = 50 germs M = 100 realisations, minimal area S = 10 pixels, maximum radius Rmax = 30 pixels). For presentation of images without zoom, the contours are dilated by a SE 3 × 3 pixels.

Figure 14 .

 14 Figure 14. First line) (a) Classification in 3 classes κ and (b) markers coming from the classification κ on image "Port de Bouc". (Second line) (d, e) deterministic approach sg mrk ( χ 2 (f λ ), κ) on (c) the chisquared metric based gradient χ 2 (f λ ). (Third line) (g, h) Stochastic approach sg mrk (mpdf κ-b (f λ ), κ) on (f ) the marginal pdf mpdf κ-b (f λ ). The parameters are Q = 7 classes, N = 50 germs M = 100 realisations, minimal area S = 10 pixels, maximum radius Rmax = 30 pixels). For presentation of images without zoom, the contours are dilated by a SE 3 × 3 pixels.

Figure 15 .

 15 Figure 15. General framework of WS controlled by a classification, with regionalized random-balls germs.The part outlined with a dotted line is not presented in this paper.

Acknowledgments

The authors are grateful to Prof. Guy Flouzat (Laboratoire de Télédétection à Haute Résolution, LTHR/ ERT 43 / UPS, Université Paul Sabatier, Toulouse 3, France) for providing the PLEIADES satellite simulated images, obtained in the framework of ORFEO program (Centre National d'Etudes Spatiales, the French space agency).

The authors would like to dedicate this paper to the memory of Prof. Guy Flouzat.