Spatial/spectral endmember extraction by multidimensional morphological operations

IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 40, NO. 9, SEPTEMBER 2002 2025 Spatial/Spectral Endmember Extraction by Multidimensional Morphological Operations Antonio Plaza, Pablo Martínez, Rosa Pérez, and Javier Plaza Abstract—Spectral mixture analysis provides an efficient mechanism for the interpretation and classification of remotely sensed multidimensional imagery. It aims to identify a set of reference signatures (also known as endmembers) that can be used to model the reflectance spectrum at each pixel of the original image. Thus, the modeling is carried out as a linear combination of a finite number of ground components. Although spectral mixture models have proved to be appropriate for the purpose of large hyperspectral dataset subpixel analysis, few methods are available in the literature for the extraction of appropriate endmembers in spectral unmixing. Most approaches have been designed from a spectroscopic viewpoint and, thus, tend to neglect the existing spatial correlation between pixels. This paper presents a new automated method that performs unsupervised pixel purity determination and endmember extraction from multidimensional datasets; this is achieved by using both spatial and spectral information in a combined manner. The method is based on mathematical morphology, a classic image processing technique that can be applied to the spectral domain while being able to keep its spatial characteristics. The proposed methodology is evaluated through a specifically designed framework that uses both simulated and real hyperspectral data. near future, the use of these types of sensors on satellite platforms will produce a nearly continual stream of high-dimensional data, and this expected high data volume will require fast and efficient means for storage, transmission, and analysis [9], [10]. Linear spectral unmixing [11]–[13] is one of the most important approaches for the analysis and classification of multi/hyperspectral datasets. This approach involves two steps: to find spectrally unique signatures of pure ground components (usually referred to as endmembers [14]) and to express individual pixels as a linear combination of endmembers [15]. Let be a spectrum of values obtained at the sensor for a certain pixel in a multispectral image. This with spatial coordinates spectrum can be considered as an -dimensional ( -D) vector (where is the number of spectral bands) and may be modeled in terms of a linear combination of several endmember vectors , according to the equations and constraints [16] (1) Index Terms—Automated endmember extraction, mathematical morphology, morphological eccentricity index, multidimensional analysis, spatial/spectral integration, spectral mixture model. (2) I. INTRODUCTION I MAGING spectroscopy, also known as hyperspectral remote sensing, allows a sensor on a moving platform to gather reflected radiation from a ground target so that a detector system consisting of CCD devices can record a great deal (typically tens) of spectral channels simultaneously. With such detail, the ability to detect and identify individual materials or classes is expected to be greatly improved [1]. During the last several years, a great deal of new airborne and spaceborne sensors have been improved for hyperspectral remote sensing applications. A chief hyperspectral sensor is the NASA/JPL Airborne Visible and Infra-Red Imaging Spectrometer (AVIRIS) system, which currently covers the wavelength region from 0.38–2.50 m [2]. Other sensors that offer hyperspectral capabilities have been available as research instruments for the last ten years, including DAIS/ROSIS [3], CASI-2 [4], AISA [5], MIVIS [6], HyMap [7], or HYPERION [8]. In the Manuscript received January 11, 2002; revised May 22, 2002. This work was supported in part by the Spanish Government under Grant TIC-2000-0739-C04-03 and by Junta de Extremadura (local Government) Grant 2PR01A085. The authors are with the Neural Networks and Signal Processing Group (GRNPS), Computer Science Department, University of Extremadura, Escuela Politécnica, 10.071 Cáceres, Spain (e-mail: aplaza@unex.es; pablomar@unex.es; rosapere@unex.es). Digital Object Identifier 10.1109/TGRS.2002.802494 is the number of endmembers needed to accurately where model the original spectrum, and is a scalar value representing . the fractional coverage of endmember vector in pixel The ideal case is that the coefficients in the linear combination are nonnegative and sum to 1, being, therefore, interpretable as cover fractions or abundances [13]. A set of endmembers and cover fractions is called a mixture model [16]. One of the new perspectives opened by this approach, together with the improved spectral resolution of sensors, is the possibility of subpixel analysis of scenes, which aims to quantify the abundance of different materials in a single pixel. Like principal components, endmembers provide a basis set in whose terms the rest of the data can be described [17], [18]. However, unlike principal components, the endmembers are expected to provide a more “physical” description of the data with no orthogonalization restrictions. For instance, a simple mixture model based on three endmembers has the geometrical interpretation of a triangle whose vertices are the endmembers. Cover fractions are determined by the position of spectra within the triangle and can be considered relative coordinates in a new reference system given by the endmembers, as Fig. 1 shows. Although the development of effective devices with new and improved technical capabilities stands out, the evolution in the design of algorithms for data processing is not as positive. Many of the techniques currently under way for hyperspectral data 0196-2892/02$17.00 © 2002 IEEE 2026 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 40, NO. 9, SEPTEMBER 2002 Fig. 1. Scatterplot of two-dimensional spectral data illustrating the physical interpretation of a mixture model based on endmembers. analysis and endmember extraction have their origins in classic spectroscopy and, thus, focus exclusively on the spectral nature of the data. Available analysis procedures do not usually take into account information related to the spatial context, which can be useful to improve the obtained classification [10]. In this paper, we introduce a novel approach to endmember extraction and multidimensional data analysis: the integration of spatial data with spectral information. Our research is structured as follows. Section II provides a detailed overview of existing endmember extraction techniques. Section III describes fundamental principles in the proposed approach. Section IV illustrates our method with examples and simulation results. Section V presents both an example of experimental data and a comparison of our method with other well-known endmember extraction approaches. Finally, Section VI includes some concluding statements and comments on plausible future research. II. BACKGROUND ON ENDMEMBER EXTRACTION TECHNIQUES A number of algorithms based on the notion of spectral mixture modeling have been proposed over the past decade to accomplish the complex task of finding appropriate endmembers for spectral unmixing in multi/hyperspectral data. One of the most successful approaches has been the pixel purity index (PPI) [15], which is based on the geometry of convex sets [19]. PPI considers spectral pixels as vectors in an -D space. A dimensionality reduction is first applied to the original data cube by using the minimum noise fraction [15]. The algorithm proceeds by generating a large number of random -D vectors, also called “skewers” [20], through the dataset. Every data point is projected onto each skewer, along which position is pointed out. The data points that correspond to extrema in the direction of a skewer are identified and placed on a list. As more skewers are generated the list grows, and the number of times a given pixel is placed on this list is also tallied. The pixels with the highest tallies are considered the purest ones, since a pixel count provides a “pixel purity index.” It is important to emphasize that the PPI algorithm does not identify a final list of endmembers. PPI was conceived not as a solution, but as a guide; in fact, the author proposed comparing the pure pixels with target spectra from a library and successively projecting the data toward lower dimensional spaces while endmembers were identified [15]. There are several interactive software tools oriented to perform this task, but the supervised nature of such tools leads to the fact that the analysis of a scene usually requires the intervention of a highly trained analyst. Then, interactive solutions may not be effective in situations where large scenes must be analyzed quickly as a matter of routine. In addition, randomness in the selection of skewers has been identified by some authors as a shortcoming of the algorithm [20]. The original implementation of PPI proposes the use of unitary vectors as skewers in random directions of the -D space. This implementation may be improved by careful selection of existing vectors to skew the dataset. Intelligent selection of skewers may result in a more efficient behavior of the algorithm. Some tools based on variations of PPI concepts have been proposed. The manual endmember selection tool (MEST) [21] is a computer display for interactive derivation of spectral endmembers. The first step in the selection process is the determination of the number of mixture components via a principal component analysis (PCA). If the PCA finds that most variance in orthogonal directions (referred the data is accounted for by to as eigenvectors), then the number of endmembers is fixed to , and MEST displays the spectral data mean corrected and projected onto the -D space, which is determined by the first eigenvectors. MEST provides the user with a means of exspectra that ploration in the mixing space to search for are acceptable as the spectral signatures of ground components; likewise, these spectra contain the spectral data in the simplex that they span in the mixing space. Again, the supervised nature of the approach may become a major drawback. There are several recent efforts in the literature that pursue the fully automated identification of endmembers from multidimensional data cubes. The N-FINDR algorithm [22] finds the set of pixels with the largest possible volume by “inflating” a simplex within the data. The input for the algorithm is the full image cube, with no previous dimensionality reduction. A random set of vectors is initially selected. In order to refine the initial estimate of endmembers, every pixel in the image must be evaluated in terms of pixel purity likelihood or nearly pure statehood. To achieve this, the volume is calculated for every pixel in the place of each endmember. A trial volume is calculated for every pixel in each endmember position by replacing that endmember and finding the volume. If the replacement results in a volume increase, the pixel replaces the endmember. This procedure is repeated until there are no more replacements of endmembers. Once the pure pixels are found, their spectra can be used to unmix the original image.This produces a set of images that show the abundance of each endmember. While the endmember determination step of N-FINDR has been optimized and can be executed rapidly, the computational performance of the algorithm depends on the accuracy of the initial random selection of endmembers. The ORASIS (optical real-time adaptative spectral identification system) algorithm [23] uses a process, called exemplar selection, to trim the dataset. This procedure rejects redundant spectra by calculating the angle between spectral vectors. Any vector that is not separated by a certain threshold angle is removed from the data. The procedure then finds a basis set PLAZA et al.: SPATIAL/SPECTRAL ENDMEMBER EXTRACTION of much lower dimension than the original data by a modified Gram–Schmidt process. The exemplar spectra are then projected onto this basis subspace, and a simplex is found through a minimum volume transform. Even though ORASIS is tuned for rapid execution, there is no dimensionality reduction as part of this algorithm. The whole process is very sensitive to the threshold angle parameter, related to sensor noise. This parameter must be selected cautiously in order to avoid the loss of important endmembers. The iterative error analysis (IEA) algorithm [24] is also executed directly on the data with neither the dimensionality reduction nor the data thinning. An initial vector (usually the mean spectrum of the data) is chosen to start the process. A constrained unmixing is then performed, and an error image is formed. The average score of vectors with higher degrees of error (distance from the initial vector) is assumed to be the first endmember. Another constrained unmixing is then performed, and the error image is formed. The average score of vectors with greater errors (distance from the first endmember) is assumed to be the second endmember. This process is continued until a predetermined number of endmembers is found, and fractional abundance maps are found by the final constrained unmixing. There are no published computer time estimates for the IEA algorithm, but the repeated constrained unmixings involved should take significant time. The use of multiple endmembers at a pixel level has also been a main topic of interest to account for endmember variability. Multiple endmember spectral mixture analysis (MESMA) [25] is a spectral unmixing approach in which many possible mixture models are analyzed simultaneously in order to produce the best fit. The weighting coefficients (fractions) for each model and each pixel are determined so that the linear combination of the endmember spectra produces the lowest margin of RMS error when compared to the apparent surface reflectance for the pixel. Weighting coefficients are constrained to lie between zero and one, and a valid fit is restricted to a maximum preset RMS error. This approach requires an extensive library of field, laboratory, and image spectra, where each plausible component is represented at least once. Since the need for spectral libraries turns out to be a demanding constraint (intensive ground work is required in their development), this fact may be the major limitation of this approach. Another attempt to incorporate endmember variability into spectral mixture analysis [16] is based on the representation of a landscape component type, not with one endmember spectrum but with a set or bundle of spectra, where each spectrum is feasible as an instance of the component (e.g., in the case of a tree component, each spectrum could be the spectral reflectance of a tree canopy). A pixel may have an infinite number of mixture models with endmembers selected from the bundles. Authors in this approach [16] use a simulated annealing algorithm to derive bundles of endmembers. This task is carried out by constructing a simplex from a partition of facets in a convex hull of a point data cloud. Bundle unmixing finds the maximum and minimum endmember fractions in those models, i.e., the range of endmember fraction values possible for a pixel. Finally, self-organizing neural network algorithms have been applied also to the unsupervised classification of multi/hy- 2027 TABLE I A COMPARISON OF SEVERAL ENDMEMBER EXTRACTION APPROACHES perspectral datasets through endmember extraction [26]. Kohonen’s self-organizing maps are suitable to perform competitive endmember learning. Hence, clusters are formed in the spectral space. The topology of the network and the training algorithm are key parameters that must be fixed with caution in order to obtain satisfactory results. Table I shows a comparison between the above-mentioned approaches, emphasizing some of their key operational features. As can be seen in Table I, most endmember extraction methods rely on spectral properties of the scene alone, neglecting the information related to the spatial correlation of pixels when deriving endmembers. It is well known that many multispectral remote sensing datasets have the unique feature of being representations of collections of images, where a particular pixel is characterized by spectral and spatial properties. Spectral information deals with pixel value in an independent manner from the neighboring pixel values, whereas spatial information concerns the relationship between the pixel value and those of its neighbors. Despite the inherent spatial–spectral duality that resides in a multidimensional data cube, most available processing techniques tend to focus exclusively on the spectral domain, handling the data not as an image but as an unordered listing of spectral measurements. These gauges can be shuffled arbitrarily without affecting analysis [10]. The incorporation of spatial information results from the notion of image processing, which is not yet fully exploited in the field of multidimensional data analysis, due to the likely fact that most analysis techniques originate from classic spectroscopy. The importance of analyzing spatial and spectral patterns simultaneously has been identified as a desired goal by many scientists who are devoted to multidimensional data analysis [27]–[29]. This type of processing can be approached from various points of view. In [30], a possibility is discussed for the refinement of results obtained by spectral techniques through a second step based on spatial context. The techniques proposed by these authors do not focus on endmember extraction, but yield some insightful feedback on plausible spatial and spectral response combinations. One of the most interesting approaches proposed is contextual classification, which accounts for the tendency of certain ground cover classes to occur more frequently in some contexts than in others. Once an initial clustering based on spectral information alone has been performed, a postclassification technique is applied. This 2028 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 40, NO. 9, SEPTEMBER 2002 method consists of two main parts: the definition of a pixel neighborhood (surrounding each pixel of the scene) and the performance of a local operation so that the pixel may be changed into the label mostly represented in the window that defines the neighborhood. This simple operation entirely separates spatial information from spectral information, and thus the two types of information are not treated simultaneously; quite the opposite, spectral information is first extracted, and a spatial context is then imposed. In this paper, we present the possibility of obtaining a method that integrates both spatial and spectral responses in a simultaneous manner. Mathematical morphology [31], [32] is a classic nonlinear spatial processing technique that provides a remarkable framework to achieve the desired integration. Morphology was originally defined for binary images and has been extended to the grayscale [33] and color [34] image cases, but it has been seldom used to process multi/hyperspectral imagery [35]. We propose a new application of mathematical morphology focused on the automated extraction of endmembers from multidimensional data. In addition, we list key advantages in the use of morphology to perform this task. A first major point is that endmember extraction is basically a nonlinear task. Furthermore, morphology allows for the introduction of a local-to-global approach in the search for endmembers by using spatial kernels (structuring elements in the morphology jargon). These items define local neighborhoods around each pixel. This type of processing can be applied to the search of convexities or extremities within a data cloud of spectral points and to the use of the spatial correlation between the data. Endmembers can be identified by following an iterative process where pixels in close proximity in the spatial domain compete against each other in terms of their convexity or spectral purity. As a result, they allow for the determination of a local representative in a neighborhood that will be compared to other locally selected pixels. The adoption of this approach leads to the incorporation of spatial information into the endmember determination procedure. As a final main step, morphological operations are implemented by replacing a pixel with a neighbor that satisfies a certain condition. In binary/grayscale morphology, the condition is usually related to the digital value of the pixel, and the two basic morphological operations (dilation and erosion) are, respectively, based on the replacement of a pixel by the neighbor with the maximum and minimum value. Since an endmember is defined as a spectrally pure pixel that describes several mixed pixels in the scene, extended morphological operations can obviously contribute to locating suitable pixels that replace others in the scene according to some desired particularity of the pixel, e.g., its spectral purity. Therefore, we propose a mathematical framework to extend mathematical morphology to the multidimensional domain, which results in the definition of a set of spatial/spectral operations that can be used to extract reference spectral signatures. Our objective in this research is to analyze the viability of extended morphological operations to analyze multidimensional remote sensing data. It is important to emphasize that the combination of spatial and spectral information leads to a new interpretation of the endmember concept, related to the capacity of a particular signature in order to describe various mixed pixels coherently in both spectral and spatial terms. III. METHODS This section is organized as follows. Section III-A provides an overview of classic mathematical morphology. Section III-B focuses on the proposed framework to extend classic morphological operations to a multidimensional domain. Some simple examples that illustrate the behavior of basic extended morphological operations are included. Finally, Section III-C describes the proposed approach for endmember extraction: the general algorithm, parameters, and implementation options are also discussed. A. Classic Mathematical Morphology Mathematical morphology is the science of shape and structure, based on set-theoretical, topological, and geometrical concepts [31], [32]. Over the past few years, mathematical morphology has grown as a powerful method for image processing and has become a solid mathematical theory. It has been successfully applied to various disciplines, such as mineralogy and medical diagnostics, finding increasing applications in digital image processing, computer vision, and pattern recognition [36]. The two basic operations of mathematical morphology are dilation and erosion. These operations were originally defined for binary images [31], but have been extended to the grayscale image case [33]. In grayscale morphology, images are handled be a gray-level function as continuous valued sets. Let representing an image, and let represent another set that comprises a finite number of pixels. This set is usually referred to as a kernel or structuring element in morphology terminology. by structuring element Erosion and dilation of image are respectively written as follows: Inf (3) Sup (4) denotes the weight associated with the different where elements of the kernel. The expressions for grayscale dilation and erosion bear a marked similarity to the convolution integral often encountered in digital image processing, with sums and differences replacing multiplication and minimum and maximum replacing summation. The main computational task in dealing with grayscale morphological operations is the query for local maxima or minima in a local search area around each image pixel. This area is determined by the size and shape of structuring element . In order to simplify the algorithms developed, we refer only to convex and plain structuring elements: a special class of struc[37]. turing elements that results when This is not a general requirement of morphological operations, and other types of kernels can be used in the future developments of the method. For illustrative purposes, let be a plain 3 3 structuring element. If a dilation operation using is applied over a grayscale image, then the local effect of the operation is the selection of the brightest pixel in a 3 3-pixel search area around the target pixel. The constraints imposed on PLAZA et al.: SPATIAL/SPECTRAL ENDMEMBER EXTRACTION 2029 the kernel definition causes all pixels lying within the kernel to be handled equally, i.e., no weight is associated with the pixels according to their position along the kernel, and the pixel with maximum digital value is selected. The previous operation is repeated with all the pixels of the image, leading to a new image (with the same dimensions as the original) where lighter zones are developed concerning the size and shape of the structuring element. In contrast, grayscale erosion has the global effect of shrinking the lighter zones of the image. The effect obtained by these operations is illustrated in Fig. 2, where a simple grayscale synthetic image is eroded and dilated using a 3 3 structuring element. By applying these basic operations, one can construct more advanced tools [31], [32]. For example, a classical morphological edge detector may be obtained by subtracting the eroded image from the dilated image; opening (erosion followed by dilation) is used to detect peaks or local maxima; closing (dilation followed by erosion) is applied to detect valleys or local minima, etc. Fig. 2. Local and global effect of applying grayscale dilation and erosion operations to a simple grayscale image. B. Extending Mathematical Morphology to Multispectral Images In grayscale morphology, the calculation of maximum and minimum gray values relies on the definition of a simple ordering relation given by the digital value of pixels. In the case of multispectral datasets, this natural assumption is not straightforward, since there is no natural means for the total ordering of multivariate pixels. Hence, the greatest challenge in the task of extending morphology operations to multispectral images is the definition of a vector ordering relation that allows for the determination of the maximum and minimum elements in any family of -D vectors [34]. In order to obtain the desired performance, we have tested two possibilities that follow a similar strategy. First, a lattice structure is imposed onto the spectral space by the definition of a distance measure, and, second, morphological operators are defined by extension. These possibilities are described in the next section. Instead, we finish this section by including a simple example of the effect of multispectral dilation and erosion as compared to grayscale counterparts. 1) Cumulative Distance for Pixel Ordering in Terms of Spectral Purity: In order to impose an ordering relation in terms of spectral purity in a set of multispectral pixels lying within a structuring element or kernel , we propose defining a metric that calculates the cumulative distance between one particular of the kernel, where denotes a vector of pixel dimensions at spatial coordinates and every other pixel . The metric is denoted by and is defined as follows: (5) is a pointwise linear distance measure between two where -D vectors, i.e., the spectral angle distance. This cumulative distance can order the vectors of a kernel in terms of their spectral purity. The computational time from calculating may be expressed in terms of the number of dot calculations as Fig. 3. Calculation of the most spectrally singular (a) and most spectrally mixed (b) pixels within a kernel by using cumulative distance . K D , where is the number of elements in the kernel. The dot operations, and distance between two vectors involves is calculated between every vector and all the others in the kernel, resulting in a cubic complexity of in asymptotical notation. Multispectral dilation and erosion can be defined after the expressions for grayscale dilation and erosion in (3) and (4), respectively, as follows: (6) (7) and , respectively, denote the pixel where vector that maximizes and minimizes cumulative distance . is a set of The ordering relation is illustrated in Fig. 3. If spatial neighboring pixels defined by a plain kernel, then is the most spectrally singular pixel [Fig. 3(a)], and is the most highly mixed element [Fig. 3(b)]. 2) Distance to the Centroid of the Data Cloud: The computational cost of calculating the cumulative distance is cubic, resulting in a very high complexity if the kernel has a large size. We have checked another possibility that relies on the calculation of the distance between every spectral point in the kernel and the centroid of the kernel data cloud, defined as (8) 2030 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 40, NO. 9, SEPTEMBER 2002 where is the number of elements in . Then, a local eccentricity measure for a particular -D point of the kernel can be obtained by calculating the distance in relation to the centroid as (9) is a pointwise linear distance function (spectral where may angle). In this case, computational time to obtain , i.e., the distance between two be expressed as dot operations and is calculated between vectors involves every vector of the kernel and the centroid. This results in a (in asymptotical quadratic computational complexity defines an ordering relation in the data notation). Distance cloud that allows for the identification of the maximum as that element that is most “eccentrical” or distant from the centroid [Fig. 4(a)], while the minimum is the closest element to the centroid [Fig. 4(b)]. Multispectral dilation and erosion can be defined respectively as Fig. 4. Calculation of the most spectrally singular (a) and most spectrally mixed (b) pixels within a kernel by using the distance to the centroid of the data cloud ( ). D K (10) (11) and , respectively, denote the pixel where vector that maximizes and minimizes the distance to the centroid . In this paper, only the two distances ( and ) are considered. We have demonstrated that the selection of the maximum and minimum through both methods is equivalent in the case of a linear distance measure as the spectral angle. Yet, this is not necessarily so. For example, if we take a nonlinear distance measure as the spectral information divergence (SID), defined in [38], the sum of the distances to each pixel in the kernel would not be equivalent, in general terms, to the distance to the centroid. Our research only focuses on a simple linear distance in order to simplify algorithm design. A comparative effort including nonlinear distances as well as distances that use second-order statistics [39], [40] deserves future attention. 3) Example of Multispectral Erosion and Dilation: The result of applying an erosion/dilation operation to a multispectral dataset is a new data cube, with exactly the same dimensions as the original, where each pixel has been replaced by the maximum/minimum of the kernel neighborhood considered. For illustrative purposes, let us consider a simple synthetic hyperspectral image with pure and mixed spectral signatures and a plain 3 3 kernel. The kernel is moved, pixel by pixel, through the image, defining a local spatial search area around each target pixel. The maximum (dilation) or minimum (erosion) pixel in the search area is calculated according to the schema described above, yielding the results depicted in Fig. 5. As can be examined in Fig. 5, dilation has the effect of expanding zones with pure spectral signatures according to the spatial properties of the structuring element, while erosion shrinks zones with pure signatures. Both operations use spatial and spectral information in a combined fashion. From this simple example, it is straightforward that extended morphological operations can be used to identify which pixels of the Fig. 5. Local and global effect of applying spatial/spectral dilation and erosion operations to a simple hyperspectral image. image scene are pure and which ones are formed by a mixture of components. C. Automated Morphological Endmember Extraction As explained above, the physical idea of an endmember is that of a spectrally pure signature (generally corresponding to a pure ground component) that can be used to describe several nonpure pixels through a mixture model based on linear combinations of endmembers [11]. The notion behind this idea is the determination of a spectral purity index at a kernel level. Two different possibilities have been explored to determine the purity index. Thus, we focus on advantages and disadvantages for both possibilities. This section ends with a detailed description of our spatial/spectral algorithm, which we have named automated morphological endmember extraction (AMEE). 1) “Voting” Approach to Endmember Selection: Multispectral dilation can extract the most particular pixel in a certain kernel neighborhood by the definition of a purity distance measure between pixels. Two distance measures have been described in Section III-B for this purpose (see Figs. 3 and 4). The identification of the purest pixel by the dilation operation is done as follows. As a first step, one particular pixel of the neighborhood is selected (candidate pixel). This pixel receives a score from all its neighbors (voters). The score is the spectral angle distance between the candidate and the neighboring pixel. Once all the neighbors have voted, the individual scores are summed, and the accumulated value is PLAZA et al.: SPATIAL/SPECTRAL ENDMEMBER EXTRACTION 2031 assigned to the candidate pixel. The operation is repeated with every pixel of the kernel, so that all the pixels assume the role of candidates once. In the end, the pixel with the highest accumulated score (the winner in the vote) is selected as the neighborhood endmember, and a counter associated with the pixel is increased in one unit. The major limitation of this approach, used by other existing endmember extraction methodologies, is the fact that such a voting system does not provide detailed information about the degree of purity in a certain pixel. For instance, the PPI algorithm uses a very similar schema. Every pixel of the scene is projected to random unitary vectors, and the most extreme pixel in the direction of the unitary vector increases its purity counter in one unit. A pixel is labeled as “spectrally pure” if it is repeatedly selected as extreme in several random projections. Nevertheless, an important question like “how extreme is the pixel in the data cloud” is neglected. The “eccentricity” of a pixel, i.e., its distance from the least spectrally pure or most highly mixed element, may be an important issue in the determination of the pixel purity degree. The incorporation of this information may result in a much more complete evaluation of the capacity of pure pixels to describe other mixed pixels and, thus, improve the whole endmember selection process. 2) “Evaluation” Approach to Endmember Selection: The voting system was characterized by a single multispectral dilation operation. In order to account for the “eccentricity” of the maximum pixel, we propose the combination of the output provided by dilation and erosion operations. While dilation selects the purest pixel (maximum), erosion selects the most highly mixed pixel (minimum) in a certain spatial neighborhood. The distance between the maximum and the minimum (a widely used operation in classic mathematical morphology) provides an eccentricity value for the selected pixel that allows for an evaluation of the quality of the pixel in terms of its spectral purity. In order to incorporate this idea, a new quality measure must be introduced. We define the morphological eccentricity index (MEI) as the distance between the maximum element (obtained by a dilation) and minimum element (obtained by an erosion) in the kernel. Fig. 6 illustrates how the MEI index is calculated. Following be the spatial coordinates of the maximum (6) and (7), let at the kernel neighborhood of a certain pixel with pixel spatial coordinates . Let be the minimum pixel of the same kernel neighborhood. Then, the MEI index is given is a pointwise distance by the following expression, where function (in this paper, the spectral angle distance). (12) The MEI index is associated with a certain pixel in order to determine its capacity to describe other pixels in the kernel. It is important to emphasize that the MEI value is assigned to the maximum pixel in the kernel neighborhood. An illustrative example of the improvements produced by the evaluation system based on the MEI index is shown in Fig. 7. The maximum element within the kernel is selected in Fig. 7(a) (voting system), but no statement can be made about its spectral purity. Fig. 6. Calculation of the MEI index by combining multispectral dilation and erosion. In Fig. 7(b) (evaluation system), the maximum element is not only selected but also evaluated using its distance to the minimum, which is particularly important when all the pixels in the kernel neighborhood have similar spectral purity. The proposed evaluation approach presents some limitations (illustrated in Fig. 8) related to the kernel size under consideration. Determining the quality of a local endmember by the MEI index is only appropriate if the kernel size is sufficiently large to include sample pixels with different degrees of purity. If the kernel used is not large enough, then there is no sufficient information to determine wether a pixel is a reasonable candidate to be used as an endmember. The main reason for that is the unavailability of different samples to perform the comparison. There may exist situations where a low MEI value is assigned to a pure pixel [as depicted in Fig. 8(a)]. In order to solve this problem, we propose the application of an iterative process with different kernel sizes. All the pixels of the scene are evaluated in each iteration, but increasingly larger neighborhoods around each pixel are considered in subsequent iterations in order to obtain detailed spatial/spectral information at a pixel level. The use of increasingly larger kernel sizes may lead to another issue. Since only one pixel per kernel neighborhood is selected, there may be situations where important endmembers are lost [as depicted in Fig. 8(b)]. This shortcoming is alleviated by the fact that the kernel is moved through all the pixels of the image, but may still produce undesired effects. Therefore, it is straightforward that the MEI index is sensitive to kernel size and shape, which must be adjusted so that a satisfactory result of the algorithm is obtained. The spatial/spectral nature of the proposed method introduces, in fact, a dependence between the structuring element and the algorithm output that is also present in classic mathematical morphology. This dependence allows the tuning of the method in order to improve the identification of 2032 Fig. 7. IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 40, NO. 9, SEPTEMBER 2002 Evaluation of endmembers improves the results obtained by voting of endmembers through the definition of a morphological eccentricity index. Fig. 8. Limitations of the evaluation approach. endmembers associated with targets with a certain size, shape, and spatial distribution in the scene. As a final concern, a minor limitation of the current MEI definition is that it tends to exaggerate the purity of a pixel in certain cases: those in which the kernel is located over an area where all the pixels are relatively pure but correspond to different endmembers. We have proved by experimenting with real data that this circumstance does not significantly affect the performance of the algorithm due to the fact that all the pixels are good candidates to be selected as endmembers. With these issues in mind, we describe the proposed endmember extraction algorithm. 3) Algorithm Description: A general block diagram of the AMEE algorithm is shown in Fig. 9. The input to the method is the full data cube, with no previous dimensionality reduction or preprocessing. Parameter refers to the number of iterations that the algorithm performs. Parameter and , respectively, denote the minimum and maximum kernel sizes that will be considered in the iterative process. These parameters are interrelated, as the number of iterations depends on the minimum and maximum kernel sizes under scrutiny. is considered. Then, the First, the minimum kernel size kernel is moved through all the pixels of the image. The spectrally purest pixel and the spectrally most highly mixed pixel are obtained at each kernel neighborhood by multispectral dilation and erosion operations. An MEI value is associated with the purest pixel by comparing the result of the dilation to the result of erosion. The previously described operation is repeated Fig. 9. Schematic representation of the proposed automated morphological endmember extraction method through a block diagram. by using structuring elements of progressively increased size, and the algorithm performs as many iterations as needed until is achieved. The associated MEI the maximum kernel size value of selected pixels at subsequent iterations is updated by means of newly obtained values, as a larger spatial context is considered, and a final MEI image is generated after iterations [41], [42]. This approach ensures a complete spatial/spectral description of the image and provides an efficient tool to integrate both types of information simultaneously. However, if the size of the maximum kernel used is large, which can become a need in order to obtain a correct description (see Fig. 7), this results in a considerable increase in morphological operation computational cost, mainly because of the high number of pixels involved in the calculation of the distances. In this case, the kernel size increases exponentially, which results in an in asymptotical notation). We exponential complexity ( have tested some variations of the previous implementation in order to reduce the computational complexity of the proposed approach, and a constrained implementation through cascade dilation operations and nonoverlapping kernel sizes [43] proved to be considerably faster than the previous approach. The constraints imposed result in a computational complexity ( in asymptotical notation). PLAZA et al.: SPATIAL/SPECTRAL ENDMEMBER EXTRACTION 2033 Fig. 10. (a) Reflectance spectra used to generate simulated scenes. (b) Radiance spectra used to generate simulated scenes. (c) Simulated Scene_1. (d) Simulated Scene_2. Automated endmember selection is performed from the MEI image obtained in the competitive endmember selection step by a threshold value . There are many available techniques for automated thresholding of grayscale images [44]. We have found appropriate results by using the multilevel Otsu method, an approach based on the minimization of the weighted sum of group variances [45]. The final selection is refined by a region-growing procedure that incorporates neighboring pixels that are sufficiently similar (according to a certain threshold angle parameter ) to the regions obtained after thresholding. Mean spectra are obtained for the resulting regions after the region-growing process. As a consequence, a final set of endmembers is obtained. This approach has been developed by following the ideas presented in [16]. Traditional methodologies such as the spectral angle mapper (SAM) or linear spectral unmixing (LSU) can be used for the purpose of mapping the obtained endmembers over the original image and obtaining a final classification result. IV. EVALUATION WITH SIMULATED DATA In this section we evaluate the performance of the proposed algorithm by using simulated data. For the sake of simplicity, we will only consider scenes with two and three different components. Section IV-A describes how synthetic data are gener- ated. Both radiance and reflectance data are used in order to analyze the effect of these two options in our endmember extraction algorithm. In Section IV-B, a preliminary experiment is conducted to study the accuracy of AMEE algorithms to classify pixels as “pure” or “not pure” and the influence of the kernel size parameter in this particular task. Finally, Section IV-C presents a set of experiments that allow us to test the accuracy of AMEE in the task of estimating the abundances of constituents. Scenes with subtle and abrupt changes in the abundance of materials are used in the simulations. A. Generation of Simulated Data A set of four spectral signatures directly extracted from a 1997 AVIRIS image scene of the Jasper Ridge, CA region [46] will be used in computer simulations to illustrate the proposed algorithm. This dataset, which has extensively been studied in the literature, is available in radiance and reflectance units. This fact allows us to study the performance of our algorithm in both cases. The selected signatures [see Fig. 10(a) and (b)] correspond to soil, vegetation, buildings, and roads. Two simple 60 60-pixel hyperspectral scenes have been generated using the above-mentioned spectra. The first scene, called Scene_1 [Fig. 10(c)], was created by using artificially generated mixtures between two selected signatures (soil and vegetation) with different abundance fractions. The scene is 2034 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 40, NO. 9, SEPTEMBER 2002 B. Identification of Pure Pixels TABLE II ABUNDANCE ASSIGNMENT FOR REGIONS IN SIMULATED SCENES SCENE_1 AND SCENE_2 formed by six regions of ten-pixels width, and the abundances have been assigned to the regions so that the linear mixture between the two components is progressive and adds to one (as shown in Table II). The synthetic hyperspectral scene in Fig. 10(d) (called Scene_2) simulates homogeneous targets (a road and a building) in a soil background. The simple simulated “objects” are a 10 10-pixel building and a road ten pixels wide. This image is characterized by abrupt changes in the abundances of materials. There is no linear mixture among components in the scene. Instead, there is a situation where 1) three different pure constituents are present and 2) each pixel is characterized by having an abundance of one for one particular endmember and zero for the others (see Table II). Both Scene_1 and Scene_2 have been generated by using radiance and reflectance spectra. We will refer to them from now on as Scene_1_radiance, Scene_1_reflectance, Scene_2_radiance, and Scene_2_reflectance. Random noise was used to simulate contributions from ambient (clutter) and instrumental sources. Noise was created by using numbers with a standard normal distribution obtained from a pseudorandom number generator and added to each pixel to generate a signal-to-noise ratio (SNR) of 30:1. For the simulations, we will define the SNR for each band as the ratio of the 50% signal level to the standard deviation of the noise [47]. Thus, the simulated data are created, based on a simple linear model, by the following expression [19]: (13) is a vector containing the simulated discrete where spectrum,; is the number of endmembers used; is a scalar value representing the fractional abundance of endmember at pixel (Table II contains a description of the abundance is the noise factor. of endmembers in each region); and The signal is scaled by 50% of the SNR, which is equivalent to reducing the noise standard deviation by the inverse factor ), so that the simulated data meet the SNR definition ( [19]. The vector terms in the parentheses are multiplied element by element. One of the properties of synthetic images is simplicity in the generation of entirely reliable ground truth from data. Receiver operating characteristic (ROC) [48] is a classic approach to evaluate the performance of image classification algorithms when ground-truth information is available. A preliminary experiment has been performed by using Scene_1_Reflectance [see Fig. 10(c)] to analyze the influence of the kernel size parameter on the accuracy of our method in the task of determining whether or not a certain pixel is “pure.” The goal of ROC analysis is to estimate a curve that describes the tradeoff between sensitivity (true positive rate) and specificity (false alarm rate) in the tested algorithm at different decision compromises [48]. In our case, the decision is related to the threshold value that determines which pixels of an MEI image, provided by the AMEE algorithm, are labeled as “pure” and which are labeled as “mixed.” Since there are many possible thresholds, the ROC curve is constructed by taking a set of equidistant threshold values between the minimum and the maximum purity index of the MEI image. Each threshold value leads to the generation of a binary classification image (BCI), where pixels with a purity index higher than the threshold are marked as “pure,” and pixels with a purity index lower than the threshold are labeled as “not pure.” Each BCI is compared to a ground-truth purity index image (GTPI) directly obtained from the original image, where only pure pixels are marked. If a certain pixel is classified as “pure” in both the BCI and GTPI, then a true positive is checked. On other hand, if a pixel is labeled as “pure” in the BCI and as “not pure” in the GTPI, then a false positive is counted. Each point of the ROC curve is obtained by plotting the average true positives against average false positives per BCI image. Fig. 11 shows an example of the construction of an ROC curve for an MEI image obtained after applying the AMEE algorithm to Scene_1_Reflectance, using a square-shaped 11 11 kernels (note that this kernel is one pixel wider than the regions in the scene). In this particular image, pure pixels are those included in regions R1 and R6. Fig. 11(a) shows the MEI image produced by the algorithm. Fig. 11(b)–(e) shows the resulting images after thresholding the MEI image with several equidistant threshold values in the dynamic range of the image. Fig. 11(f) shows an example of ROC curve construction by using the previous images. Each point at the curve is obtained by calculating the average number of true and false positives per image. As can be pinpointed in the graph, the extreme points of an ROC curve are (0, 0) and (1, 1). These values indicate that the price of complete success in the algorithm for recognizing all existing targets is complete failure, i.e., detection of all false alarms as true positives. Then, as the performance of a method increases, the ROC curve shifts toward the upper left corner, which represents perfect performance. Then, the quality of an MEI image may be determined by calculating the area below its associated ROC curve. From detection theory definitions [48], the greater the area below the ROC curve, the better the algorithm performance. In order to analyze the impact of the kernel size parameter on the classification, we have performed the following experiment. PLAZA et al.: SPATIAL/SPECTRAL ENDMEMBER EXTRACTION 2035 2 Fig. 11. (a) MEI image obtained after applying the AMEE algorithm to Scene_1_Reflectance by using a square-shaped 11 11 kernel. (b)–(e) Resulting images after thresholding the MEI image with four equidistant threshold values. (f) Construction of an ROC curve for the MEI image (the area under the curve is also addressed). (g) Area under ROC curves associated with classification images obtained by using voting and evaluation of endmembers with different kernel sizes. The kernel size parameter is initially set to 3 3 pixels and is progressively increased to a maximum value of 21 21 pixels. For each considered size, an MEI image is obtained. An ROC curve is then constructed for each MEI image, and the area under the curve is estimated. The resulting area estimations are plotted against the correspondent values of the parameter, obtaining a curve that indicates which values of the parameter provide better results. Fig. 11(g) shows the results of this experiment after applying the AMEE algorithm to Scene_1_Reflectance. For comparative purposes, we have considered the voting and evaluation systems for endmember selection (refer to Section III for a detailed description of each approach). The results shown in Fig. 11(g) reveal that the purity determination process is sensitive to the kernel size parameter, which must be large enough to contain samples with different spectral purity. In both cases (voting and evaluation), performance is high when the kernel size is 11 11 or larger (the width of spatial patterns in Scene_1 is ten pixels). In contrast, this experiment reveals that voting is much more dependent on kernel size than evaluation, which provides acceptable results for small kernel sizes. C. Estimation of Abundances In order to perform abundance estimation simulations, we have applied our algorithm to the simulated scenes described in Fig. 10. A set of endmembers is extracted for each image, and the abundance of each endmember is estimated by using fully constrained linear spectral unmixing [49]. In order to determine the accuracy of our method in this particular task, we compare the estimated abundances to the true abundances, shown in Table II. Both reflectance and radiance spectra are used in the simulations. Fig. 12 shows some results obtained after applying the AMEE algorithm to Scene_1_reflectance and Scene_1_radiance using a square-shaped 11 11 kernels. In both cases, two endmembers were extracted. The reflectance endmembers are labeled as Ref_1_1 and Ref_1_2 [Fig. 12(a)] and the radiance endmembers as Rad_1_1 and Rad_1_2 [Fig. 12(b)]. A visual comparison of the extracted endmembers to the original signatures used to generate the images [Fig. 10(a) and (b)] reveals that the extracted endmembers correspond to soil and vegetation. Even though the resulting abundance maps are noisy [Fig. 12(c)–(f)], the spatial zones corresponding to pure soil and vegetation can be easily identified. A scatter plot of true versus estimated fractional abundances is also shown for each abundance map [Fig. 12(g)–(j)], showing that the correlation is high even in the presence of noise. It should also be pointed out that the algorithm produces very similar results with radiance and reflectance data. This is not a surprising result: the transformation of radiance data into reflectance units can modify the spectral shape, but the spectral singularity of different constituents remains. Fig. 13 shows the results obtained after applying the AMEE algorithm to Scene_2_reflectance and Scene_2_radiance. Three endmembers were extracted for each scene [Fig. 13(a) and (b)]. A visual comparison of the extracted endmembers to the original signatures used to generate the images [Fig. 10(a) and (b)] reveals that the extracted endmembers correspond to soil, road, and building. By using the road and building extracted endmembers, we have created two matched filters and their correspondent output images [Fig. 13(c)–(f)]. Even though the output is noisy, we can easily separate the two targets from the background. These objects can be easily obtained from the images through simple thresholding [the 2036 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 40, NO. 9, SEPTEMBER 2002 Fig. 12. (a) Endmembers derived from Scene_1_Reflectance. (b) Endmembers derived from Scene_1_Radiance. (c), (d) Abundance maps for reflectance endmembers. (e), (f) Abundance maps for radiance endmembers. (g), (h) Scatterplot of true versus estimated abundance for reflectance endmembers. (i), (j) Scatterplot of true versus estimated abundance for radiance endmembers. automated Otsu method [45] was used to obtain the threshold value, producing the images shown in Fig. 13(g)–(j)]. Some of the above-mentioned segmented images are corrupted by “salt-and-pepper” noise due to the sensor noise simulation in the scenes. This noise can be easily eliminated, at the expense of some image degradation, by simple 3 3 median filters [50]. The experiments presented in this section reveal that the AMEE algorithm can be used for the purpose of mixture quantification and abundance determination. The algorithm is also applicable for the purpose of identifying and separating targets from a homogeneous background. The use of radiance/reflectance data does not significantly affect the performance of the method. Although the method seems to work well in theory development, further experiments with real data are desirable in order to fully validate the previous remarks. V. EVALUATION WITH EXPERIMENTAL DATA In this section, a comparative analysis of AMEE to other well-known endmember extraction algorithms is presented. The comparison has been carried out with real AVIRIS hyperspectral data from the Cuprite mining district in Nevada. This scene is well understood mineralogically [51] and has reliable ground truth in several forms, including laboratory spectra [52] and abundance maps [53]. This fact has made this region a standard test site for comparison of classification algorithms. The scene that we have selected for this work was acquired in 1995 and consists of 50 spectral bands in the SWIR-2 region of the spectrum, where the differences among minerals are clear. In order to perform this experiment, we have selected several laboratory spectra from the U.S. Geological Survey (USGS) Digital Spectral Library [52] to be used as ground-truth references in the comparison. These signatures, acquired PLAZA et al.: SPATIAL/SPECTRAL ENDMEMBER EXTRACTION 2037 Fig. 13. (a) Endmembers derived from Scene_2_reflectance. (b) Endmembers derived from Scene_2_radiance. (c), (d) Matched filter output for building and road (Scene_2_reflectance). (e), (f) Matched filter output for building and road (Scene_2_radiance). (g), (h) Otsu segmentation of (c) and (d). (i), (j) Otsu segmentation of (e) and (f). with ground spectrometers on the Cuprite site, correspond to minerals that have been identified in the 1995 AVIRIS scene. It is important to emphasize that there are other minerals present in the scene, as well as intimate mixtures between several minerals, but we have focused our study on four particular minerals (alunite, buddingtonite, calcite, and kaolinite) that can be found prominently and in pure form in the area. Two of the most well known algorithms for endmember extraction are compared to the AMEE algorithm in this section: PPI and N-FINDR. As described in Section II, PPI requires close human supervision during the process of determining the endmembers, while N-FINDR is fully automated. PPI and N-FINDR endmembers were extracted and provided by the developers of each method. The parameters used by the authors in order to produce the above-mentioned results are unknown to us. The kernel size parameter used by the AMEE algorithm was 15 15. This parameter was determined after analyzing the width in pixels of patterns of interest in the abundance planes published in [53]. The compared spectra have been labeled as follows. The USGS spectra were labeled with their correspondent names in the library. AMEE, N-FINDR, and PPI endmembers were labeled with the initial letter of the method plus a number that indicates the order of identification of the signatures. Each method has its own mechanisms to extract the endmembers, and the order in which the endmembers are identified is not important. Thus, AMEE endmembers are labeled as A1…A9, N-FINDR endmembers as N1…N10, and PPI endmembers as P1…P11. It should be taken into account that the number of endmembers is different for each method and that “false” endmember spectra are produced due to shade, noise, sensor artifacts, etc. The term “false” should not be misunderstood, since these endmembers are actually needed to accurately unmix the scene. In order to determine the best matching endmember for each USGS reference signature, the confusion matrix method [54] is used. Tables III–V show, respectively, the confusion matrix of spectral angle distances between USGS ground-truth 2038 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 40, NO. 9, SEPTEMBER 2002 TABLE III CONFUSION MATRIX OF SPECTRAL ANGLE DISTANCES BETWEEN AMEE ENDMEMBERS AND USGS REFERENCE SPECTRA TABLE IV CONFUSION MATRIX OF SPECTRAL ANGLE DISTANCES BETWEEN N-FINDR ENDMEMBERS AND USGS REFERENCE SPECTRA TABLE V CONFUSION MATRIX OF SPECTRAL ANGLE DISTANCES BETWEEN PPI ENDMEMBERS AND USGS REFERENCE SPECTRA signatures and AMEE, N-FINDR, and PPI endmembers. We have selected the spectral angle distance for the comparison due to existing scale and ilumination variations between reference and extracted endmembers, which are mainly due to atmospheric transmission effects that are not present in spectra acquired by ground spectrometers. The endmember spectra for AMEE, N-FINDR, and PPI are shown in Fig. 14 for the following minerals: 1) 2) 3) 4) alunite; buddingtonite; calcite; kaolinite. Extracted endmembers have the same reflectance units as the input data, since each endmember is actually the spectrum (or an average) of real image pixels. For illustrative purposes, ground USGS spectra for the same minerals are also plotted. USGS spectra are expressed in scaled reflectance. It can be noted that, although there are differences in the scales, spectral shape is preserved in most cases. A visual comparison of the buddingtonite USGS laboratory spectrum and derived endmembers [Fig. 14(b)] reveals that atmospheric transmission effects are relevant in the case of this mineral. Although the three algorithms find endmembers by using very different techniques, the output of the three algorithms operating on the same dataset is very similar, as expected from the results shown in the confusion matrices. For a quantitative analysis of the performance of each approach, the plots in Fig. 14 also address the spectral angle distance between each extracted endmember and its correspondent USGS reference spectrum. AMEE endmembers present high spectral similarity to ground truth for the alunite, buddingtonite, and kaolinite minerals. N-FINDR provides the best result for the calcite mineral. The combination of spectral and spatial context information in the pure pixel identification step of AMEE is the main reason for the encouraging results obtained by the algorithm in this example. VI. CONCLUSION AND FUTURE LINES Spectral unmixing and object classification are two important goals of remote imaging spectroscopy. The idea of using endmembers derived from the data for classification and unmixing has been considered before, but few methods have exploited the existing correlation between neighboring pixels. The integration of spatial and spectral information has been identified as a desired goal by the scientific community dedicated to multidimensional data analysis and classification. We have described a novel approach based on mathematical morphology to perform unsupervised endmember extraction from multi/hyperspectral data. The method uses both spectral and spatial information in a simultaneous manner. An evaluation of the method by using simulated and real data has been presented and discussed. PLAZA et al.: SPATIAL/SPECTRAL ENDMEMBER EXTRACTION 2039 Fig. 14. Algorithm scene derived and USGS laboratory reflectance spectra for (a) alunite, (b) buddingtonite, (c) calcite, and (d) kaolinite. Extracted endmembers have the same reflectance units as the input data, while USGS spectra are expressed in scaled reflectance. The spectral angle distance between each extracted endmember and the correspondent USGS reference spectrum is addressed. Results with simulated data reveal that the method can accurately model the spatial distribution of spectral patterns in the scene by extended morphological operations that apply plain spatial kernels. The spatial properties (size and shape) of the kernel have a strong influence on the final result, a fact that is consistent with classic mathematical morphology theory. The behavior of the algorithm also depends on the relationship between the spatial properties of the kernels used and the distribution of spectral patterns in the scene. This fact allows for the tuning of our method for a wide range of applications, from target detection to global classification and spectral unmixing of scenes. The use of reflectance/radiance data does not have a significant impact on the output of the algorithm. Results with experimental data show that the proposed method produces results that are comparable to those found by working with other widely accepted methodologies. The proposed method is accurate in the task of identifying endmembers from complicated scenes as the famous AVIRIS Cuprite dataset, which has become a standard test site for the comparison of algorithms due to the availability of quality ground-truth measurements. As with any new approach, there are some unresolved issues that may present challenges over time. In this sense, future lines will cover some relevant developments that were not included in the present study. An evaluation of different distance measures (both linear and nonlinear) to be used in the extension of morphological operations is a key topic, deserving future research. In addition, the structuring elements developed in this paper are limited to plain square-shaped kernels. The use of kernels with no such limitations is of great interest in order to explore, in greater detail, the existing spatial/spectral interrelation between the kernel used and patterns in the scene. The number of endmembers extracted per kernel neighborhood is also an interesting issue. In the present research, some limitations of the proposed method have been identified as a consequence of the fact that only one pixel per kernel neighborhood is selected. Alternative definitions of the local MEI index may overcome existing limitations. Finally, efficient hardware implementations based on field-programmable logic arrays and systolic arrays are being currently tested at our laboratory in order to provide the methodology with real-time capabilities. ACKNOWLEDGMENT The authors acknowledge the suggestions and comments of J. Pinzón, R. O. Green, and J. A. Gualtieri that helped to improve 2040 IEEE TRANSACTIONS ON GEOSCIENCE AND REMOTE SENSING, VOL. 40, NO. 9, SEPTEMBER 2002 the quality of this paper. Indications made by P. Gomez-Vilda about efficient implementation in the proposed method are also gratefully acknowledged. In addition, the authors would like to express their gratitude to M. Winter and E. Winter for providing results of the N-FINDR algorithm for the Cuprite dataset. Finally, we also wish to state our appreciation to A. Curado Fuentes, from the Department of English at University of Extremadura, who reviewed the submitted version of this paper. REFERENCES [1] N. M. Short et al.. The remote sensing tutorial. NASA/GSFC. [Online]. Available: http://rst.gsfc.nasa.gov. [2] R. O. Green et al., “Imaging spectroscopy and the airborne visible/infrared imaging spectrometer (AVIRIS),” Remote Sens. Environ., vol. 65, pp. 227–248, 1998. [3] A. Müller, A. Hausold, and P. Strobl, “HySens—DAIS/ROSIS imaging spectrometers at DLR,” in Proc. SPIE Image and Signal Processing for Remote Sensing VII, Toulouse, France, 2001. [4] I. Keller and J. Fischer, “Details and improvements of the calibration of the compact airborne spectrographic imager (CASI),” in First EARSEL Workshop on Imaging Spectroscopy, Paris, France, 1998, pp. 81–88. [5] C. Mao, M. Seal, and G. Heitschmidt, “Airborne hyperspectral image aquisition with digital CCD video camera,” in 16th Biennial Workshop on Videography and Color Photography in Resource Assessment, Weslaco, TX, 1997, pp. 129–140. [6] R. Bianchi, R. M. Cavalli, L. Fiumi, C. M. Marino, and S. Pignatti, “CNR LARA project, Italy: Airborne laboratory for environmental research,” in Summaries of the V JPL Airborne Earth Science Workshop, Pasadena, CA, 2001. [7] P. Launeau, J. Girardeau, D. Despan, C. Sotin, J. M. Tubia, B. Diez-Garretas, and A. Asensi, “Study of serpentine and specific vegetation in the Ronda Peridotite (Spain) from 1991 AVIRIS to 2001 hymap images,” in Summaries of the XI JPL Airborne Earth Science Workshop, Pasadena, CA, 2001. [8] F. A. Kruse, J. W. Boardman, and J. F. Huntington, “Progress report: Geologic validation of E0-1 hyperion using AVIRIS,” in Summaries of the XI JPL Airborne Earth Science Workshop, Pasadena, CA, 2001. [9] S. Tadjudin and D. Landgrebe, “Classification of high dimensional data with limited training samples,” Ph.D. dissertation, School of Elect. Eng. Comput. Sci., Purdue Univ., Lafayette, IN, 1998. [10] V. Madhok and D. Landgrebe, “Spectral–spatial analysis of remote sensing data: An image model and a procedural design,” Ph.D. dissertation, School of Elect. Eng. Comput. Sci., Purdue Univ., Lafayette, IN, 1998. [11] J. J. Settle, “On the relationship between spectral unmixing and subspace projection,” IEEE Trans. Geosci. Remote Sensing, vol. 34, pp. 1045–1046, July 1996. [12] Y. H. Hu, H. B. Lee, and F. L. Scarpace, “Optimal linear spectral unmixing,” IEEE Trans. Geosci. Remote Sensing, vol. 37, pp. 639–644, Jan. 1999. [13] M. Petrou and P. G. Foschi, “Confidence in linear spectral unmixing of single pixels,” IEEE Trans. Geosci. Remote Sensing, vol. 37, pp. 624–626, Jan. 1999. [14] F. A. Kruse, “Spectral identification of image endmembers determined from AVIRIS data,” in Summaries of the VII JPL Airborne Earth Science Workshop, Pasadena, CA, 1998. [15] J. W. Boardman, F. A. Kruse, and R. O. Green, “Mapping target signatures via partial unmixing of AVIRIS data,” in Summaries of the V JPL Airborne Earth Science Workshop, Pasadena, CA, 1995. [16] C. A. Bateson, G. P. Asner, and C. A. Wessman, “Endmember bundles: A new approach to incorporating endmember variability into spectral mixture analysis,” IEEE Trans. Geosci. Remote Sensing, vol. 38, pp. 1083–1094, Mar. 2000. [17] J. W. Boardman and F. A. Kruse, “Automated spectral analysis: A geological example using AVIRIS data, Northern Grapevine Mountains, Nevada,” in Proc. 10th Thematic Conference, Geologic Remote Sensing, San Antonio, TX, 1994. [18] F. A. Kruse, J. F. Huntington, and R. O. Green, “Results from the 1995 AVIRIS geology group shoot,” in Proc. 2nd Int. Airborne Remote Sensing Conf. Exhib., 1996. [19] A. Ifarraguerri and C.-I. Chang, “Multispectral and hyperspectral image analysis with convex cones,” IEEE Trans. Geosci. Remote Sensing, vol. 37, pp. 756–770, Mar. 1999. [20] J. Theiler, D. Lavenier, N. Harvey, S. Perkins, and J. Szymanski, “Using blocks of skewers for faster computation of pixel purity index,” in SPIE Int. Conf. Optical Science and Technology, San Diego, CA, 2000. [21] C. A. Bateson and B. Curtiss, “A tool for manual endmember selection and spectral unmixing,” in Summaries of the V JPL Airborne Earth Science Workshop, Pasadena, CA, 1993. [22] M. E. Winter, “N-FINDR: An algorithm for fast autonomous spectral end-member determination in hyperspectral data,” in Proc. SPIE Imaging Spectrometry V, 1999, pp. 266–275. [23] J. Bowles, P. J. Palmadesso, J. A. Antoniades, M. M. Baumback, and L. J. Rickard, “Use of filter vectors in hyperspectral data analysis,” Proc. SPIE Infrared Spaceborne Remote Sensing III, pp. 148–157, 1995. [24] K. Staenz, T. Szeredi, and J. Schwarz, “ISDAS—A system for processing/analyzing hyperspectral data,” Can. J. Remote Sens., vol. 24, pp. 99–113, 1998. [25] D. Roberts, M. Gardener, J. Regelbrugge, D. Pedreros, and S. Ustin, “Mapping the distribution of wildfire fuels using AVIRIS in the Santa Monica Mountains,” in Summaries of the VIII JPL Airborne Earth Science Workshop, Pasadena, CA, 1998. [26] P. Martínez, J. A. Gualtieri, P. L. Aguilar, R. M. Pérez, M. Linaje, J. C. Preciado, and A. Plaza, “Hyperspectral image classification using a selforganizing map,” in Summaries of the XI JPL Airborne Earth Science Workshop, Pasadena, CA, 2001. [27] J. Pinzón, S. Ustin, and J. F. Pierce, “Robust feature extraction for hyperspectral imagery using both spatial and spectral redundancies,” in Summaries of the VII JPL Airborne Earth Science Workshop, Pasadena, CA, 1998. [28] L. Hudgins and C. Hines, “Spatial–spectral morphological operators for hyperspectral region-growing,” in Int. Symp. Spectral Sensing Research (ISSSR), Quebec City, Canada, 2001. [29] B. Peppin, P. L. Hauff, D. C. Peters, E. C. Prosh, and G. Borstad, “Different correction-to-reflectance methods and their impacts on mineral classification using SFSI and low- and high-altitude AVIRIS,” in Summaries of the XI JPL Airborne Earth Science Workshop, Pasadena, CA, 2001. [30] L. O. Jimenez-Rodriguez and J. Rivera-Medina, “Integration of spatial and spectral information in unsupervised classification for multispectral and hyperspectral data,” Proc. SPIE Image and Signal Processing for Remote Sensing V, pp. 24–33, 1999. [31] J. Serra, Image Analysis and Mathematical Morphology. London: Academic Press, 1982. , Image Analysis and Mathematical Morphology. London: Aca[32] demic Press, 1993, vol. 1. [33] S. R. Sternberg, “Greyscale morphology,” Comput. Vision, Graphics Image Processing, vol. 35, pp. 283–305, 1999. [34] P. Lambert and J. Chanussot, “Extending mathematical morphology to color image processing,” in Proc. CGIP 2000, Saint-Etienne, France, 2000. [35] P. Soille, “Morphological partitioning of multispectral images,” J. Electron. Imaging, vol. 5, pp. 252–265, July 1996. [36] R. M. Haralick and L. G. Shapiro, Computer and Robot Vision. Reading, MA: Addison-Wesley, 1993. [37] M. Comer and E. Delp, “Morphological operations for color image processing,” J. Electron. Imaging, vol. 8, pp. 279–289, July 1999. [38] C.-I. Chang, “An information theoretic-based approach to spectral variability, similarity and discriminability for hyperspectral image analysis,” IEEE Trans. Inform. Theory, vol. 46, Aug. 2000. [39] L. O. Jimenez-Rodriguez and D. Landgrebe, “Hyperspectral data analysis and supervised feature reduction via projection pursuit,” IEEE Trans. Geosci. Remote Sensing, vol. 37, pp. 2653–2667, Nov. 1999. , “Supervised classification in high-dimensional space: Geomet[40] rical, statistical and asymptotical properties of multivariate data,” IEEE Trans. Syst., Man, Cybern. A, vol. 28, Feb. 1998. [41] A. Plaza, P. Martínez, J. A. Gualtieri, and R. M. Pérez, “Spatial/spectral endmember extraction from AVIRIS hyperspectral data using mathematical morphology,” in Summaries of the XI JPL Airborne Earth Science Workshop, Pasadena, CA, 2001. [42] P. L. Aguilar, A. Plaza, R. M. Pérez, and P. Martinez, “Morphological endmember identification and its systolic array design,” in Neural Networks and Systolic Array Design, Singapore: World Scientific Publishers, 2002, ch. 3. [43] A. Plaza, P. Martinez, J. A. Gualtieri, and R. M. Perez, “Automated identification of endmembers from hyperspectral data using mathematical morphology,” in Proc. SPIE Image and Signal Processing for Remote Sensing VII, Toulouse, France, 2001. [44] P. K. Sahoo, S. Soltani, and A. K. C. Wong, “A survey of thresholding techniques,” Comput. Vision, Graphics Image Processing, vol. 41, pp. 233–260, 1988. PLAZA et al.: SPATIAL/SPECTRAL ENDMEMBER EXTRACTION [45] N. Otsu, “A threshold selection method from gray-level histograms,” IEEE Trans. Syst., Man, Cybern., vol. SMC-9, Jan. 1979. [46] California Institute of Technology. AVIRIS (Airborne Visible/Infrared Imaging Spectrometer) homepage. [Online]. Available: http://aviris.jpl.nasa.gov. [47] J. C. Harsanyi and C.-I. Chang, “Hyperspectral image classification and dimensionality reduction: An orthogonal subspace projection approach,” IEEE Trans. Geosci. Remote Sensing, vol. 32, pp. 779–785, July 1994. [48] N. A. Macmillan and C. D. Creelman, Detection Theory: A User’s Guide. Cambridge, U.K.: Cambridge University Press, 1991. [49] D. Heinz and C. I Chang, “Fully constrained least squares linear mixture analysis for material quantification in hyperspectral imagery,” IEEE Trans. Geosci. Remote Sensing, vol. 39, pp. 529–545, Mar. 2001. [50] R. C. Gonzalez and R. E. Woods, Digital Image Processing. Reading, MA: Addison-Wesley, 1992. [51] G. Swayze, “The hydrothermal and structural history of the cuprite mining district, Southwestern Nevada: An integrated geological and geophysical approach,” Ph.D. dissertation, Univ. Colorado, Boulder, CO, 1997. [52] R. N. Clark, G. A. Swayze, A. Gallagher, T. V. V. King, and W. M. Calvin, “The U.S. Geological Survey Digital Spectral Library: Version 1: 0.2 to 3.0 m,”, U.S. Geological Survey Open File Rep. 93–592, 1993. [53] R. N. Clark and G. A. Swayze, “Evolution in imaging spectroscopy analysis and sensor signal-to-noise: An examination of how far we have come,” in Summaries of the VI JPL Airborne Earth Science Workshop, Pasadena, CA, 2001. [54] R. G. Congalton, “Considerations and techniques for assessing the accuracy of remotely sensed data,” in Proc. IGARSS, vol. 3, 1989, pp. 1847–1850. Antonio Plaza received the M.S. degree in computer science from the University of Extremadura, Caceres, Spain, in 1997, where he is currently pursuing the Ph.D. degree. He is currently an Assistant Professor at the Computer Science Department, University of Extremadura, since 1999, where he was a Research Associate, 1997 to 1999. He was also a Visiting Researcher at Applied Information Sciences Branch, Goddard Space Flight Center, Greenbelt, MD, and the AVIRIS Data Facility, Jet Propulsion Laboratory, Pasadena, CA. His main research interests include computer vision, image processing, pattern recognition, and the development and efficient implementation of algorithms for multi/hyperspectral image analysis. 2041 Pablo Martínez received the M.S. and Ph.D. degrees in physics from the University of Granada, Andalucia, Spain, in 1980 and 1992 respectively. He is currently the Head Scientist of the Neural Networks and Signal Processing Group (GRNPS) at the University of Extremadura, Caceres, Spain, where he has also been a Professor of computer science at University of Extremadura since 1985. He was Visiting Researcher at the Applied Information Sciences Branch, Goddard Space Flight Center, Greenbelt MD in 2000. His main research interests include remote sensing, digital image analysis, and neural networks. Rosa Pérez received the M.S. degree in mathematics from the University of Extremadura, Caceres, Spain, in 1985, and the Ph.D. degree in computer science from the Polytechnic University of Madrid, Madrid, Spain, in 1995. Since 1985, she has been a Professor of computer science at University of Extremadura, Caceres, Spain. Her main research interests include neural networks, systolic arrays, field-programmable logic array design, and signal processing. Javier Plaza is currently pursuing the M.S. degree in computer science from the University of Extremadura, Caceres, Spain. His current research work is focused on imaging spectrometry. He is mainly involved in quantitative applications using real and simulated remote sensing data. Other major research interests include telecommunications and image processing.