A survey of methods incorporating spatial information in image classification and spectral unmixing

ABSTRACT

Over the past decade, the incorporation of spatial information has drawn increasing attention in multispectral and hyperspectral data analysis. In particular, the property of spatial autocorrelation among pixels has shown great potential for improving understanding of remotely sensed imagery. In this paper, we provide a comprehensive review of the state-of-the-art techniques in incorporating spatial information in image classification and spectral unmixing. For image classification, spatial information is accounted for in the stages of pre-classification, sample selection, classifiers, post-classification, and accuracy assessment. With regards to spectral unmixing, spatial information is discussed in the context of endmember extraction, selection of endmember combinations, and abundance estimation. Finally, a perspective on future research directions for advancing spatial-spectral methods is offered.

1. Introduction

Image classification and spectral unmixing are two predominant methods used for multispectral and hyperspectral data analysis. Along with the advancement of both satellite remote sensing and spatial science, spatial information has attracted significant interest within the remote-sensing community. Numerous studies have been conducted with the aim of effectively incorporating spatial autocorrelation in both image classification and spectral unmixing. However, compared to spectral information, there are few summaries of recent efforts to investigate how spatial information is manifested and their effectiveness in remote-sensing analysis. In this paper, we provide a comprehensive review of state-of-the-art techniques used to incorporate spatial information in image classification and spectral unmixing.

2. Incorporating spatial information in image classification

In this section, we summarize the state-of-the-art methods that incorporate spatial information in the image classification process according to the following five categories: (1) pre-classification, (2) sample selection, (3) classifiers, (4) post-classification, and (5) accuracy assessment.

2.1. Pre-classification

2.1.1. Image texture

2.1.1.1. Grey-level co-occurrence matrix (GLCM)

The grey-level co-occurrence matrix (GLCM), introduced by Haralick, Shanmugam, and Dinstein (1973), is probably the most commonly adopted texture extraction method for remote-sensing images. A GLCM records the frequency of pairs of pixel values, where the pairs are oriented in a particular direction (0°, 45°, 90°, 135°), that occur in an image. Various texture measures, such as the mean, variance, homogeneity, contrast, dissimilarity, entropy, second moment, and correlation, can then be calculated from this matrix. The GLCM (and thus the texture measures) are usually calculated within a moving window to generate the texture images. While many studies have explored different GLCM texture measures and the role of different window sizes to better characterize image texture, very few have attempted to improve the statistical method itself. Clausi and Zhao (2002) used a hash table and linked list (grey-level co-occurrence hybrid structure (GLCHS)) to store the co-occurring probabilities, resulting in a significant reduction in the computation time for calculating the co-occurrence texture features. Later, the GLCHS algorithm was again improved using a grey-level co-occurrence hybrid histogram (GLCHH) for sum and difference histograms (Clausi and Zhao 2003). By integrating GLCHS and GLCIA, the resulting grey-level co-occurrence integrated algorithm (GLCIA) requires only 0.04–16% of the computational time of the standard GLCM method to determine co-occurrence probability texture features (837–850). More recently, Huang, Liu, and Zhang (2014a) applied clustering and sparse representation to multispectral images to generate codes or quantized levels based on the multi-channel texture features to be extracted. As a result, the traditional single-band GLCM method was improved so that multi-/hyperspectral imagery features can be extracted. Similarly, Tsai and Lai (2013) were able to extract texture features from image cubes rather than 2-D neighbourhoods by extending the original co-occurrence matrix to a 3-D tensor field.

2.1.1.2. Geostatistical approach

Geostatistical methods have been widely applied to characterize the spatial relationships among pixel values of a particular remote-sensing image (i.e. texture) (Atkinson and Lewis 2000). The variogram is a summary of how semivariance (a measure of autocorrelation) varies as a function of distance between image pixels. There are, in general, two approaches to extracting texture features using geostatistical methods: 1) using the actual values of the sample variograms directly (Berberoglu et al. 2007; Carr and De Miranda 1998; Maillard 2003); or 2) using coefficients extracted from modelled variograms (Balaguer et al. 2010; Herzfeld and Higginson 1996). Traditionally, both approaches extract texture features individually from a single band and thus information only from the selected band or bands is used. However, these methods can be extended using multivariate variograms in multispectral dimensions (Li, Cheng, and Guo 2009). In a recent study, it was found that cross-variograms derived from different bands of the same image or same bands from different images can also improve classification results when seasonal/temporal information is important (Rodriguez-Galiano et al. 2012). In addition, compared to conventional per-pixel methods, sample variograms (and thus the texture features) can also be extracted from image objects/fields (Balaguer et al. 2010). This object-oriented approach was found to be more effective for characterizing spatial variations in certain landscapes (e.g. Mediterranean land cover) where between-field spatial variations confound the within-field spatial variations that is of interest (Lloyd et al. 2004).

2.1.1.3. Fractal approach

Fractal geometry is used to describe the irregular and highly complex shapes of natural objects such as coastlines and landscapes, which the traditional Euclidean geometry cannot analyse. The fractal dimension (often denoted D), is a non-integer value that measures the irregularity of complex objects. A perfectly flat surface has a fractal dimension of 2.0 and the value of D will increase and approach 3.0 as the surface begins to fill a volume. Therefore, features in a remotely sensed image usually have a fractal dimension of any number between 2.0 and 3.0. Based on the assumption that pixels having the same fractal dimension share the same physical properties, fractal dimensions have been used as a texture measure in remote-sensing image classification (De Jong and Burrough 1995; Emerson, Lam, and Quattrochi 2005; Sun 2006). However, there is a growing consensus that methods used to derive fractal texture images suffer from many limitations and drawbacks, and fractal dimension alone is not sufficient to provide complete texture information (Myint 2003; Parrinello and Vaughan 2002; Sun et al. 2006). One example is that the moving window used for local fractal dimension calculation must be large enough to perform a regression while small window sizes are desired for better classification accuracies. Silvetti and Delrieux (2013) proposed a solution to this problem by computing the local fractal dimension based on quadratic self-correlation and showed that this improved method worked well even with very small windows. Similarly, combining other texture measures or using multifractal dimensions were also proposed to overcome the limitation of using one single fractal dimension. Pant, Singh, and Srivastava (2010) used both fractal dimension and spatial autocorrelation statistics (Moran’s I) and the combined texture image yielded better classification results. Improved classification accuracies were also achieved by adopting multi-fractal textural descriptors, since multi-fractal measures consider the whole spectrum of dimensions and thus are able to provide a higher level of textural information (Parrinello and Vaughan 2006; Wawrzaszek, Krupinski, and Aleksandrowicz. 2013).

2.1.1.4. Filter-based approach

Although filter-based approaches have been extensively applied for image interpretation in computer science, it was not until recently that these methods were introduced to the field of remote sensing for texture feature extraction purposes (Jain and Healey 1998). Gabor filters and wavelet transformation, among many other image-filtering methods, are the two methods most frequently adopted in the remote-sensing literature. Both methods are well known for their ability to characterize image texture at multiple scales by incorporating either different spatial frequencies (Gabor filters) or multiple levels of decomposition (wavelet transformation) (Angelo and Haertel 2003; Myint, Lam, and Tyler 2004; Zhu and Yang 1998). However, textures were extracted only from single-band grey-level images in early applications. The introduction of the opponent features allowed Gabor filters to characterize texture information across multiple spectral bands (Jain and Healey 1998; Shi and Healey 2003). Similarly, Huang, Zhang, and Li (2008) extracted both spectral and spatial feature sets from remote-sensing images by applying separate wavelet transformation on multispectral bands and PCA band 1. More recently, the spectral dimension was added to traditional spatial dimensions in both filters. It was found that, when image filtering was conducted in the spectral-spatial domain, both the 3-D Gabor filters (Bau, Sarkar, and Healey 2010; Shen and Jia 2011; Shen et al. 2013) and 3-D wavelet transformation (Guo, Huang, and Zhang 2014; Huang and Zhang 2012) performed better than traditional 2-D approaches, as the new methods were capable of capturing correlations in space, wavelength, and joint spatial/spectral correlation.

2.1.2. Morphological image analysis

Another set of methods that can incorporate spatial information at the pre-classification stage is morphological image analysis, or mathematical morphology (MM) (Soille and Pesaresi 2002). MM can extract or suppress desired image objects by using set operators (such as union, intersection, and complementation) with a set of carefully selected (for both shape and size) structuring elements (SEs). The two most widely used morphological operators are opening and closing, which are combinations of the fundamental erosion and dilation operators. Aside from deleting objects smaller than the SEs, opening and closing also modify or deform objects. To avoid this problem and preserve the shape of objects, reconstruction is often used. Opening or closing by reconstruction preserves the shapes of image objects (pixels and their neighbouring pixels) by using connected operators (such as geodesic operators). However, Bellens et al. (2008) showed that due to the fact that many connected simple objects in remote-sensing images (such as roads and buildings) are treated as a single object by the reconstruction process, this often leads to unexpected and undesirable results. This problem of ‘over-reconstruction’ can be solved by partial reconstruction, which limits the extent of reconstruction (Liao et al. 2012).

When applying MM in remote-sensing applications, it is also important to include a range of SEs with different sizes to characterize different structures present in the image. This can be addressed by using morphological profiles (MPs), which are defined using successive opening/closing operations with an SE of increasing size (Pesaresi and Benediktsson 2001). The resulting MP consists of both an opening profile and a closing profile with a range of different SE sizes, and thus is a powerful tool for extracting spatial information from remotely sensed images. Nevertheless, MP cannot be used in multispectral dimensions. The most commonly adopted method of MP construction for multispectral images is the Extended Morphological Profile (EMP) (Benediktsson, Palmason, and Sveinsson 2005). An EMP is constructed by building MPs using spectral features extracted from the original image (via feature extraction methods such as principle component analysis or independent component analysis). Therefore, an EMP contains both spectral information from the extracted features and spatial information enclosed in the MPs. Classification accuracy of high spatial resolution hyperspectral data using EMPs increased significantly compared to traditional spectral classification (Fauvel et al. 2008, 2013). More recently, in order to overcome limitations such as the fixed shape of SEs, the lack of ability to describe grey-level-related information, and the high computational complexity of the MPs and EMPs (Ghamisi, Benediktsson, Cavallaro, et al. 2014), Dalla Mura et al. (2010) proposed the concept of morphological attribute profiles (APs). In contrast to morphological filters used in MPs, APs are constructed by the sequential application of morphological attribute filters (AFs). As a result, image decomposition (and thus attributes of a connected component) can be based on any type of measurement, such as scale (which is used in MPs), shape, and texture (Ghamisi, Benediktsson, and Sveinsson 2014; Ghamisi, Mura, and Benediktsson 2015). Similar to EMPs, extended APs (EAPs) have been used to handle multispectral/hyperspectral images.

2.1.3. Other spatial information extraction methods at pre-classification stage

Spatial autocorrelation or spatial statistics, such as Moran’s I and Geary’s C, have been used in several studies to characterize spatial correlation among pixels in remote-sensing images (Emerson, Lam, and Quattrochi 2005; Myint 2003; Myint and Lam 2005). These statistics are often calculated based on the spectral variability of the pixels within a predefined window (with a certain size and shape), making it difficult to extract geometric information regarding image objects. Chen, Qin, et al. (2014) proposed a directional spatial correlation (DSC), based on Moran’s I and Geary’s C, to characterize the size, shape, and spectral variability of image objects. The algorithm delineates the geometry of an object by expanding from the centre pixel along eight different directions with various lengths. Therefore, objects of different shapes and sizes can be characterized by the calculated measures. Similarly, Zhang et al. (2006) defined a pixel shape index (PSI) for each image pixel to describe shape features by using extended direction lines around a central pixel. The performance of PSI was later further improved by Huang, Zhang, and Li. (2007) using a structural feature set (SFS) extracted from the directional lines histogram (including PSI).

2.2. Sample selection

It was established many years ago that spatial autocorrelation should be taken into account when sampling reference data, which is then used as either training or validation data (Congalton 1991; Campbell 1981). Because geographical data have spatial autocorrelation, the assumption of independence of the samples may be violated if groups of contiguous pixels are selected. Therefore, the design of sample unit, sample size, and sampling scheme should account for spatial information.

2.2.1. Sample unit

The potential sampling units have been grouped into three categories: 1) a single pixel; 2) a cluster of pixels; and 3) an object (Congalton and Green 2009; Chen and Stow 2002). When the area of each sample unit increases, the positional error has diminished impact on the reported accuracy.

Single pixels are used as training samples in order to avoid the effects of spatial autocorrelation. However, it has been found that using single pixels does not guarantee a higher classification accuracy than using clusters of pixels. On the other hand, there may be the disadvantage of under-representing the spatial information for spatially heterogeneous classes (Chen and Stow 2002). When used as validation samples, single pixels have also been preferred less because of the difficulty in registering the ground-collected samples with the classification image at single-pixel accuracy (Congalton and Green 2009).

The spatial autocorrelation of the training data, using clusters of pixels, was found to decrease when the spatial resolution of the image was coarser (Chen, Stow, and Gong 2004). Small clusters of pixels have been found more effective in representing the spatial and spectral information for classes with high spatial heterogeneity (Chen and Stow 2002). The validation data can also be collected at a unit of clustered pixels (e.g. 3 pixels × 3 pixels) instead of a single pixel to mitigate the bias in the estimate of the accuracy caused by the spatial misregistration between the classification result and the validation data (Foody 2008). The pixels within each cluster or block are all used as training samples. A sensitivity analysis by Mountrakis and Xi (2013) confirmed that these spatially grouped validation samples can achieve similar classification precision to that of the spatially random samples.

It has been claimed that for object-based classification, objects, but not pixels, should be used as sampling units (Congalton and Green 2009; Stehman and Wickham 2011). Nevertheless, Radoux et al. (2011) after reviewing 20 OBIA-related papers found that only 40% of them used objects as sampling units, and that a majority of research are still using pixels as sampling units.

Objects as sample units have advantages for object-based classification. First, the sample unit corresponds to the mapping units. Second, potential error due to spatial misregistration is greatly reduced (Congalton and Green 2009). On the other hand, when objects are used as sample units for validation data, there are some critical issues: 1) heterogeneity may exist within a sample unit (Stehman and Wickham 2011); and 2) the sample units in one validation data set may have various sizes and thus areas (Radoux and Bogaert 2014). Accuracy assessment methods aiming to resolve these issues have been developed (Radoux et al. 2011; Lizarazo 2014). It has been found that larger objects are better classified (Radoux and Bogaert 2014; Castilla et al. 2014). Therefore, it is important to have objects of different sizes in the validation sample set (1029–1037).

2.2.2. Sample size

The optimal sample size of the training data is the minimum sample size that is required to achieve a certain classification accuracy (Foody et al. 2006). (Chen and Stow 2002), after conducting experiments using images of different spatial resolution levels (2 m, 4 m, 8 m, 12 m, and 16 m), found a positive relation between the required number of training samples and spatial heterogeneity.

The optimal sample size of the validation data is the minimum sample size that is required to estimate the classification accuracy with a certain precision (Congalton and Green 2009; Baraldi, Boschetti, and Humber 2014). Traditional algorithms for theoretically determining the optimal sample size are mainly based on binomial or multinomial distributions, aiming at reaching a pre-specified class-specific accuracy and confidence interval (Congalton 1991; Baraldi, Boschetti, and Humber 2014). Nevertheless, the optimal sample size is dependent on the classification scenario and the sampling design. It has been found that the optimal sample size for object-based accuracy assessment is smaller than that of pixel-based accuracy assessment (Radoux et al. 2011; MacLean et al. 2013).

2.2.3. Sampling scheme

Although probability sampling protocols such as simple random sampling or systematic random sampling are readily available, in practice it is common to select samples from specific areas where the true label of the ground is known (Zhen et al. 2013).

Stehman (2009) suggested seven criteria for the sampling sheme for selecting reference samples. One of these is that the samples should be evenly distributed in space in order to represent the entire study area. To ensure spatial balance, a common protocol for sampling, especially when the study area is large, is to first define spatially evenly distributed subareas, and then randomly select sample units from each of the subareas (Zhao et al. 2014; Dihkan et al. 2013).

A typical approach is to separate a reference data set into training and validation samples (Mountrakis and Xi 2013). In that case, it’s important to consider the spatial dependence between the training and validation samples. If both data sets are randomly selected within the same patches, there tends to be an optimistic bias (Zhen et al. 2013). Therefore, it is recommended to separate the known patches into training patches and validation patches before random sampling.

2.3. Classifiers

2.3.1. Object-based image analysis

Object-based image analysis (OBIA) paves the way for combining spectral and spatial information, and in doing so potentially offers a more comprehensive classification approach. It seeks to create meaningful image objects by segmenting an image into groups of adjacent pixels with similar characteristics (e.g. spatial and spectral). It uses image objects rather than just individual pixels for image analysis (Benz et al. 2004; Hay and Castilla 2008). With the advent and proliferation of the availability of high-resolution imagery, the use of OBIA has increased rapidly (Blaschke 2010; Blaschke et al. 2014; Liu, Guo, and Kelly 2008). OBIA primarily encompasses two stages: image segmentation and classification. In the image segmentation process, an image is partitioned into non-overlapping homogeneous regions (i.e. image objects) based on one or more criteria of homogeneity (e.g. colour and shape in the software package eCognition) (Fu and Mui 1981; Tarabalka et al. 2012; Blaschke 2010). In general, within-object variance is minimized compared to between-object variance (Desclée, Bogaert, and Defourny 2006). Each segmented image object can be characterized by both spectral and spatial properties (Jensen 2000; Benz et al. 2004; Blaschke and Strobl 2001; Baatz and Arno 2000). Spectral properties of an image object are typically statistical measures of spectral values, such as the pixels’ mean, standard deviation, and spectral ratios. Spatial properties include the shape (e.g. area, length, compactness, and smoothness), neighbourhood, and topologies. Spatial properties explicitly embodied by the individual image object are crucial to OBIA and are the main reason for the marked increase in OBIA-related studies in recent years (Hay et al. 2003; Blaschke et al. 2014; Blaschke, Lang, and Hay 2008; Hay and Castilla 2006; Conchedda, Durieux, and Mayaux 2008). During the classification stage, homogeneous image objects, instead of individual pixels, are analysed based on their spectral and spatial properties using traditional parametric or machine learning classification approaches, knowledge-based approaches, or fuzzy logic approaches (Jensen 1996; Jensen 2000; Benz et al. 2004). An iterative application of segmentation and classification is usually desirable in OBIA for extracting the information of interest (Benz et al. 2004). The initial classification of segmented image objects can serve as the high-level input for subsequent segmentation (i.e. classification-based segmentation) to improve the classification accuracy.

Over the last decade, a variety of segmentation methods have been developed to partition an image into relatively homogenous image objects (Fauvel et al. 2013; Pal and Pal 1993; Neubert, Herold, and Meinel 2008). From an algorithmic perspective, these methods can be categorized into three groups: point-based, edge-based, and area-based (Pal and Pal 1993; Schiewe 2002). Depending on the features employed in the segmentation process, the methods can be categorized into the following groups: spatial-based, spectral-based, and spatial-spectral-based (Fauvel et al. 2013; Fu and Mui 1981). Spatial-based methods search for homogenous regions of spatially connected pixels with the defined criteria (e.g. greyscale texture). Representative methods in this group include region growing, split-merge methods, and watershed methods (Blaschke, Burnett, and Pekkarinen 2004; Gonzalez, Woods, and Eddins 2004). Spectral-based methods, such as thresholding and partitional clustering, group pixels into spectral clusters based on spectral similarity measures, without considering the spatial locations of these pixels. The disjointed regions within the same cluster are assigned different labels (Fu and Mui 1981; Fauvel et al. 2013; Tarabalka, Benediktsson, and Chanussot 2009; Celeux and Govaert 1992). Spatial-spectral-based methods exploit both spatial and spectral information in the segmentation process, and allow for the merging of spatially disjoint regions (Tilton 1998). One example is HSeg segmentation, which is based on hierarchical stepwise optimization (HSWO) (Beaulieu and Goldberg 1989; Tilton et al. 2012). HSeg segmentation combines region growing with unsupervised classification to merge the pre-segmented image objects in terms of spectral similarity measures iteratively, and generates a hierarchical sequence of image objects with different levels of detail (Fauvel et al. 2013; Plaza et al. 2009; Tarabalka et al. 2012; Tilton et al. 2012). Technical details of image segmentation methods are reviewed in Pal and Pal (1993), Fu and Mui (1981), and Ilea and Whelan (2011) .

One critical parameter in implementing the segmentation is scale (Blaschke 2010; Benz et al. 2004; Burnett and Blaschke 2003; Wang, Sousa, and Gong 2004). The larger the scale value, the larger the image objects grow in the segmentation (Baatz and Schäpe 1999). Though a single scale (i.e. one-level representation) might be sufficient or straightforward for representing the content in specific applications (Tiede, Lang, and Hoffmann 2008), it is often not possible to determine a single appropriate scale in advance. Besides, image objects of interest usually have their own inherent scales. Significant image objects may be generated at different scales of analysis (Hay et al. 2001, 2003; Woodcock and Strahler 1987). Multi-scale representation is thus more appropriate for scene understanding and image interpretation. However, this multi-scale representation cannot be achieved simply by changing the resolution of the imagery, since this results in the loss of much useful context information (Benz et al. 2004). OBIA conducted at multiple scales can overcome this limitation to a large degree. Hence, establishing a hierarchical multi-scale image object network is important for appropriately interpreting the phenomena and understanding relations within an image (Hay and Castilla 2008; Burnett and Blaschke 2003). This hierarchical image object network can be built with eCognition, which builds on fractal net evolution segmentation and is currently the most widely used object-based image analysis system (Benz et al. 2004). Through defining the scale range and homogeneity criteria, image objects can be created at multiple scales in eCognition wherein each successive scale contains larger segmented image objects (Baatz and Arno 2000; Baatz et al. 2000; Kim et al. 2011). To date, a significant number of research studies (e.g. 50–55% of more than 800 articles identified by Blaschke (2010)) were conducted with eCognition software. Yet the strong dependence of the OBIA research on eCognition software dilutes the influence of various alternative hierarchical segmentation approaches on object classification accuracy. A wide spectrum of hierarchical segmentation approaches (e.g. hierarchical split merge refinement (HSMR), multi-scale object-specific analysis (MOSA), and HSeg) need to be investigated and evaluated to construct a more comprehensive object-based image analysis system (Hay et al. 2001; Tilton et al. 2012; Ojala and Pietikäinen. 1999).

At each segmentation scale, the image object in the network is characterized by its spectral-spatial features and neighbours. It is also topologically linked to its sub-objects and super-objects in a strict hierarchical structure (Benz et al. 2004; Mallinis et al. 2008). The hierarchical scale dependencies among the image objects allow for understanding scaling relations and exploring substructures of general land-cover regions. Methods for upscaling image objects hierarchically have been investigated, specifically for evaluating whether segmented image objects created at various scales should be built precisely one upon another, or generated independently with the delineated outlines fitting exactly (Blaschke 2010; Hay et al. 2005, 2001; Castilla, Hay, and Ruiz-Gallardo 2008). Owing to the hierarchical relations, different procedures can be implemented in different types of super-objects (e.g. urban and agricultural areas) to more appropriately characterize the landscape structure and explain the phenomenon (Benz et al. 2004). Moreover, the hierarchical relationship between sub-objects and super-objects can be used to determine the validity of sub-object classes (e.g. pasture), given the super-object type (e.g. urban areas). As the magnitude of aggregation and abstraction, scale is a crucial aspect of image understanding and a hot research topic in OBIA (Blaschke 2010; Blaschke et al. 2014).

Geographic object-based image analysis (GEOBIA) has recently been perceived as a new and evolving paradigm in remote sensing and geographic information science (GIScience) (Blaschke et al. 2014). Compared to OBIA, GEOBIA emphasizes the geographical space of image objects and the integration of image objects with GIS applications. GEOBIA provides a spatial perspective for understanding various physical and abstract image objects (e.g. different agricultural systems in different locations), with the ultimate goal of developing spatial knowledge or geographic-based intelligence (Blaschke et al. 2014; Blaschke 2010). GEOBIA typically involves the human semantics, ontologies, and perceptions as guiding principles for describing complex scene content, and conceptualizing knowledge based on meaningful image objects within a geographical context (Lang 2008). This spatial expert knowledge incorporation seems to be an increasing trend in remote sensing and GIScience (Blaschke et al. 2014; Andres, Arvor, and Pierkot 2012; Yue et al. 2013), but might not be so obvious to computer vision, signal processing, and pattern recognition research that focuses on advancing segmentation or classification algorithms. More cooperation across disciplines needs to be encouraged to promote the development of this multidisciplinary research field.

2.3.2. Pixel-based classifiers

In recent years, spatial information has been increasingly explored in the pixel-based classification process, especially for high-resolution imagery. Most methods investigate the role of spatial information (e.g. texture) as additional feature bands without modifying the classifiers. Contextual classifiers, however, exploit the spatial information among neighbouring pixels in classifiers to improve the classification results. Two representative methods, namely contextual support vector machines and Markov random fields, are discussed below. Sparse representation, which has recently been explored as a method to incorporate spatial information, is then reviewed. Finally, as a means to quantify the spatial structure and dependence, geostatistical methods are also reviewed.

2.3.2.1. Contextual support vector machines

Over the last decade, intensive efforts have been conducted to advance the performance of classifiers in remote sensing (Camps-Valls et al. 2014; Fauvel et al. 2013; Plaza et al. 2009). Among all classifiers, support vector machines (SVMs) are particularly appealing to the remote-sensing community (Mountrakis, Im, and Ogole 2011). SVMs have generally been found not to be susceptible to the Hughes phenomenon (i.e. a degradation of classification performance with increasing numbers of input variables) and have shown remarkable performance in numerous applications, especially with limited training samples. A thorough review can be found in Mountrakis, Im, and Ogole (2011). As a non-parametric pixel-wise statistical learning technique, SVM aims to construct an optimal separating hyperplane that maximizes the margin (i.e. space) between classes based on support vectors. Support vectors are defined as the nearest vectors from each class to the hyperplane (Suykens and Vandewalle 1999; Vapnik 2013; Cortes and Vapnik 1995). When data cannot be linearly separated, SVM solves the inseparability problem with kernel functions, which project the data into a higher dimensional feature space (Bioucas-Dias et al. 2013; Cristianini and John 2000; Schölkopf and Smola 2002). The commonly used kernel functions in SVMs are Gaussian radial basis function (RBF), polynomial, linear, and sigmoid functions (Chang and Lin 2011; Schölkopf and Smola 2002). Conventionally, these kernels are employed to project the spectral vectors of pixels into the higher dimensional feature space, without explicitly considering the spatio-contextual information.

Owing to the properties of kernel functions (e.g. a scaled summation of Mercer’s kernels are valid kernels), several SVM variants have recently been developed to simultaneously exploit the spatial and spectral information (Plaza et al. 2009; Camps-Valls et al. 2014; Fauvel et al. 2013). One approach is composite kernel-based SVM, which includes a weighted linear combination of spectral and spatial kernels computed with available features (Camps-Valls et al. 2006; Fauvel, Chanussot, and Benediktsson 2012; Fauvel et al. 2013). Spatial kernels define the projection rule of pixel-based spatial features, which can be derived with statistical filtering (e.g. mean, median, and variance) applied to each pixel’s neighbourhood (e.g. morphological neighbourhood or fixed neighbourhood) (Fauvel et al. 2013; Fauvel, Chanussot, and Benediktsson 2012). Compared to stacked spatial and spectral information as feature bands in SVM, the composite kernel allows specification of the relative influence of the extracted features, and offers more flexibility for integrating diverse information in land-cover mapping (Fauvel, Chanussot, and Benediktsson 2012; Camps-Valls et al. 2014). However, spatial kernels in the composite currently only summarize basic information (e.g. median) of the spatial neighbourhood. More sophisticated techniques for describing spatial structures (e.g. shape, size, and homogeneity) of the neighbourhood are called for (Fauvel, Chanussot, and Benediktsson 2012). Recently, this composite kernel has been further extended to multiple kernel learning, which provides a kernel framework for fusing multi-source spectral and spatial information in SVM classification (Tuia et al. 2010a, 2010b).

The introduction of composite kernels opens a wide range of opportunities for exploring the use of spatial information in classifiers. Instead of extracting spatial features before kernel computation, a recent advance in SVM is the use of graph kernels to capture all possible spatial relations between pixels of a neighbourhood in the graph-based feature space (Camps-Valls, Shervashidze, and Borgwardt 2010). Through working directly with the original data, graph kernel-based SVM overcomes the limitations associated with explicit spatial feature computation and is also advantageous in exploiting spatial relations at multiple scales.

2.3.2.2. Markov random fields

Markov random fields (MRFs) are a powerful family of stochastic models for exploring the spatial dependence properties associated with images on a per-pixel basis. They generalize the concepts of a 1-D Markov chain to a 2-D framework (e.g. pixel lattices) for characterizing contextual information and exploiting the continuity among neighbouring pixels in a classifier (Moser, Serpico, and Benediktsson 2013; Plaza et al. 2009). In the Bayesian decision theory framework, the class label of each pixel is predicted through the Bayesian maximum a posteriori (MAP) rule, which optimizes the parameters in the framework by maximizing the posterior mass functions (Moser, Serpico, and Benediktsson 2013; Li 2009; Kindermann and Snell 1980). However, due to the large number of parameters to be optimized jointly, this global MAP image classification is computationally intractable.

MRFs offer a computationally feasible solution to this issue by assuming that all the pixels satisfy the Markovianity condition. This condition implies that the probability distribution of the class label of each pixel is only affected by the labels of its neighbouring pixels (Li 2009; Geman and Geman 1984). With this assumption, MRFs can significantly reduce the computational burden in optimizing the parameters in the Bayesian framework, by changing from a global model complying with the MAP rule, to a local model based on its neighbourhood system. MRF-based approaches have therefore gained increasing popularity in remote sensing over the last decade (Li 2009; Moser, Serpico, and Benediktsson 2013; Melgani and Serpico 2003; Tso and Olsen 2005; Jia and Richards 2008).

In MRFs, based on the Hammersley–Clifford theorem, the Bayesian MAP rule can be formulated as the minimization of the energy function (Geman and Geman 1984; Dubes and Jain 1989). The basic energy function can be expressed as a linear combination of two parts. One part is related to the conditional pixel-wise spectral information, and it is the same as the non-contextual Bayesian classification system. The second part (also known as Potts model or multi-level logistic model) considers the spatio-contextual information of each pixel with the Kronecker function, which favours the generation of connected regions in spatially homogeneous land-cover classes (Jackson and Landgrebe 2002; Moser, Serpico, and Benediktsson 2013; Geman and Geman 1984). To minimize the energy function in random fields, two classical iterative approaches have been adopted in the literature, namely simulated annealing (SA) and iterated conditional mode (ICM) (Geman and Geman 1984; Solberg, Taxt, and Jain 1996; Moser, Serpico, and Benediktsson 2013). However, the long computation time of SA and the local convergence issue associated with ICM limit their abilities in energy minimization. A number of global and local minimization methods, including message passing (e.g. loopy belief propagation and tree-reweighted message passing), convex relaxation, and semi-global labelling, have been proposed for tackling the optimization problem (Szeliski et al. 2008; Yedidia, Freeman, and Weiss 2003; Kolmogorov 2006; Wainwright, Jaakkola, and Willsky 2005; Schindler 2012). One approach that has recently drawn considerable attention is graph-cut (Kolmogorov and Zabih 2004; Boykov and Funka-Lea 2006; Greig, Porteous, and Seheult 1989; Boykov, Veksler, and Zabih 2001). This method minimizes the energy function through solving a maximum flow problem over a graph. Graph-cut methods have been demonstrated to converge to the global minimum for two classes, or a strong local minimum for more than two classes, with acceptable computation time (Kolmogorov and Zabih 2004).

MRFs exhibit considerable capability in incorporating multiple spatio-contextual information sources in the classification process. The information is not limited to the spatial dependence among neighbouring pixels, but also involves spatial edges, geometrical structures, multi-scale analysis, and spatial anisotropy (Moser, Serpico, and Benediktsson 2013; Li 2009; Solberg, Taxt, and Jain 1996; Moser and Serpico 2008). The energy function in MRFs can be extended to a linear combination of multiple contribution parts, with each information source being one part. This flexibility in the choice of energy function makes MRFs a powerful tool for meeting the specific requirements of complex land-cover classification systems, and taking into account the various roles of spatial information in diverse applications to improve classification accuracy (Moser, Serpico, and Benediktsson 2013; Moser and Serpico 2008, 2011). MRFs have therefore demonstrated advantages over other pixel-based classification methods. A current challenge with MRFs is to develop more novel and reliable methods, such as the integration of MRFs with SVMs, to extract accurate remotely sensed information (Zhang et al. 2011; Moser and Serpico 2013). Investigating new input information sources (e.g. new types of spatial information) also plays a crucial role in advancing this spatio-contextual classifier in land cover mapping.

As a generative model, MRFs comply with the MAP rule through modelling the joint probability of the observed data and the corresponding class labels (i.e. the product of class conditional probability (also named likelihood) and class prior probability) (Li 2009). MRFs need to explicitly model the distribution of the likelihood and assume that the observed image data are conditionally independent given the class labels (Kumar and Hebert 2003; Besag 1974). This conditional independence assumption makes the computation tractable, but meanwhile limits the contextual information to the labelled image (Kumar and Hebert 2003; Zhong and Wang 2010). To directly take into account the spatio-contextual information in both observed and labelled images, conditional random fields (CRFs) have been proposed and recently utilized in image classification (Lafferty, McCallum, and Pereira 2001; Wallach 2004; Kumar and Hebert 2003; Zhong and Wang 2010). CRF assumes that class labels, conditional on the observed data, obey the Markov property. As a discriminative model, CRF directly models the posterior probability of class labels, given the observed image data, as the Gibbs distribution. It avoids the problem of explicit modelling of the likelihood in MRFs and relaxes the corresponding conditional independence assumption. CRFs can be used to classify various remotely sensed images regardless of the distributions they follow, and are capable of incorporating the spatial contextual information in both observed and labelled images (Zhao, Zhong, and Zhang 2015; Zhong, Zhao, and Zhang 2014; Kosov, Rottensteiner, and Heipke 2013; Hoberg et al. 2015; Zhong, Lin, and Zhang 2014).

Similar to MRFs, the maximization of the posterior probability in CRFs can be formulated as the minimization of the energy function. The most commonly used CRF energy function in image classification can be written as the sum of the unary and pairwise potentials (Kumar and Hebert 2003). The unary potential is formulated to model the pixel-wise relationship between the spectral information of each pixel and its corresponding class label without considering its neighbours. The minimization of the energy that consists of only this potential is equivalent to the non-contextual Bayesian classification of the image. Generalized linear models (e.g. multinomial logistic regression) have been utilized in the unary potential to model the class posteriors given the observations (Kumar and Hebert 2003; Zhong and Wang 2010). Yet the unary potential is not restricted to a specific form. Recently, state-of-the-art discriminative classifiers, such as support vector machines, Gaussian mixture models, and random forests, have also been explored (Kosov, Rottensteiner, and Heipke 2013; Schindler 2012; Zhong, Lin, and Zhang 2014; Zhao, Zhong, and Zhang 2015). The pairwise potential is used to model the spatial contextual information of each pixel in terms of its neighbourhood in both observed and labelled images. The pairwise potential learns a data-dependent pairwise smoothness term (e.g. generalized Ising model) with an assumption of smoothness to favour the neighbourhood pixels with similar class labels in homogeneous image regions (Zhong and Wang 2010; Kumar and Hebert 2003). It takes into account the neighbouring class labels and the whole observed image data in determining the class label of each pixel. CRFs have been demonstrated to perform better than non-contextual classification methods (e.g. random forests), but may oversmooth the image with the relatively simple smoothness term. More complex and sophisticated smoothness functions thus need to be investigated in CRFs (Schindler 2012). Recently, distance boundary constraints have been explored in the pairwise potential to preserve the spatial details in the image classification result (Zhang and Jia 2012; Zhong, Lin, and Zhang 2014). Additionally, CRFs can incorporate high-order potentials, which can model more wide-ranging and complex contextual information with the high-order spatial neighbourhood (Zhong and Wang 2011). Yet these high-order potentials increase the complexity of models and make efficient inference problematic.

2.3.2.3. Sparse representation

Sparse representation has recently been investigated in classifying remotely sensed imagery. In sparse representation, each pixel in the image is assumed to be represented by a sparse linear combination of training samples from a dictionary (Chen, Nasrabadi, and Tran 2011; Zhang et al. 2014b; Song et al. 2014). The dictionary can be designed based on the mathematical model of the data, or a set of realizations of the data (Rubinstein, Bruckstein, and Elad 2010). The coefficients associated with training samples in sparse representation are the sparse vector, whose non-zero entries denote the weights of the selected training samples. The sparse vector can be recovered through sparsity-constrained optimization, which aims to minimize the approximation error between the observed and reconstructed pixel values, given the sparsity level constraint. The constraint of the sparsity level is used to set the upper bound of the number of non-zero entries in the sparse vector. In sparse representation, the spatial contextual information from the neighbouring pixels can be exploited through imposing another constraint (i.e. smoothness constraint) in the optimization process. For example, two types of smoothness constraints, namely the Laplacian constraint and the joint sparsity model, have been devised within the sparsity-constrained optimization formulation (Chen, Nasrabadi, and Tran 2011). The Laplacian constraint forces the reconstructed pixel of interest to share similar spectral characteristics with those of its four nearest neighbours, while the joint sparsity model assumes that the pixels in a small neighbourhood can be modelled simultaneously by the sparse linear combination of a few common training samples. The sparse vector recovered through optimization can then be utilized to determine class labels of pixels of interest. Recently, kernel sparse representation has been proposed to represent pixels in terms of sparse linear combination of training samples in a feature space induced by a kernel function (Chen, Nasrabadi, and Tran 2013; Liu et al. 2013). In the optimization formulation, the kernelized joint sparsity model, for example, is served as the smoothness constraint in the induced feature space to incorporate the contextual information across neighbouring pixels (Chen, Nasrabadi, and Tran 2013). Kernel sparse representation has been demonstrated to improve data separability and outperform the conventional linear sparse representation.

2.3.2.4. Geostatistical methods

Rooted in variography, geostatistical techniques offer a valuable means to quantify the spatial image structure and spatial dependence between observations (Curran and Atkinson 1998; Atkinson and Lewis 2000). Variograms are key functions in geostatistics as they are essential in characterizing the spatial variance structure, with the stationary assumption, and play a crucial part in other geostatistical techniques (e.g. kriging and stochastic simulation) (Atkinson and Lewis 2000; Van Der Meer 2012). In recent years, variogram-based texture measures from images have been increasingly used as additional feature bands to be incorporated into the classification process (Atkinson and Lewis 2000; Van Der Meer 2012). Geostatistical techniques, however, are still rarely used as classifiers or spatial variants in classifiers (Atkinson and Naser 2010). Among the few studies conducted, indicator kriging has been used to exploit the spatial variability of classes within the image and estimate class probabilities for each pixel (Van Der Meer 1994). This geostatistical method has been shown to be more accurate than conventional classification approaches (Das and Singh 2009; Park, Chi, and Kwon 2003). Additionally, geostatistically weighted K-nearest neighbour (k-NN) has been employed to incorporate the spatial variation between classes for classifying remotely sensed imagery. This classifier modifies the original k-NN with three geographical weightings and allows for the integration of the feature information from both spatial and spectral domains (Atkinson and Naser 2010; Atkinson 2004).

2.4. Post-classification

Spectral-based pixel-wise classified images usually suffer from salt-and-pepper noise, due to the lack of spatial information in the classification (Lu and Weng 2007). Spatio-contextual information can be incorporated into the post-classification process to reduce the noise effect and build more spatially consistent thematic maps (Fauvel et al. 2013). The spatio-contextual post-processing typically consists of two elements: spatial neighbourhoods and spatial regulation. Spatial neighbourhoods can be obtained through image segmentation, or simply arbitrarily defined (e.g. first- and second-order neighbourhoods) (Tarabalka, Benediktsson, and Chanussot 2009; Fauvel et al. 2013; Tarabalka et al. 2010b). Spatial regulation can then be performed to modify the label of each individual pixel based on its neighbourhood. The regulation process can help filter out spatially isolated pixels, enhance spatial patterns of classified images, and create more spatially homogenous maps. The most commonly used spatial regulation method is majority voting, which assigns the pixels to the most frequent class in its spatial neighbourhood. Majority voting has been used in several studies to refine pixel-wise classification (e.g. SVMs) results in segmented spatial neighbourhoods (Tarabalka, Chanussot, and Benediktsson 2010; Tarabalka, Benediktsson, and Chanussot 2009; Tarabalka et al. 2010a; Fauvel et al. 2013). Additionally, Gaussian smoothing, bilateral filtering, edge-aware filtering, and anisotropic diffusion have been explored as spatial regulation methods to locally smooth classification results, through taking into account the labels, probabilities, and intensity values of the pixels within the spatial neighbourhood (Schindler 2012; Huang et al. 2014b).

Spatial reclassification kernel (SPARK) has also been employed to regularize classified images (Barnsley and Barr 1996). SPARK refines the class of each centred pixel in defined kernels (or fixed window sizes) based on the frequency and spatial arrangements of adjacent land-cover classes. The kernel size for each individual classified pixel in SPARK is further optimized with a similarity index to take into account the spatial variation of different land-cover types (Van Der Kwast et al. 2011). SPARK exhibits potentials for distinguishing complex spatial assemblages and subtle differences of land-use types in urban areas (Barnsley and Barr 1996). Moreover, relearning methods (i.e. relearning-Hist and relearning-PCM) have recently been proposed to adaptively learn the frequency and spatial arrangement of class labels in the neighbourhood (Huang et al. 2014b). Relearning methods have been found efficient in correcting classification errors and enhancing class separability.

MRF-based spatial regulation methods have recently been explored to incorporate spatio-contextual information in the post-classification process. Through minimizing the global energy function, MRF-based spatial regulation methods can iteratively refine the initial spectral pixel-wise classification results to obtain the MAP estimate (Tarabalka et al. 2010b). Unlike the original spectral image used in MRF-based classifiers (see Section 2.3.2), the observed data in MRF-based spatial regulation represent the probabilistic pixel-wise classification result (Tarabalka et al. 2010b; Huang et al. 2014b). The MRF-based spatial regulation method is advantageous, with its flexibility to exploit multiple spatio-contextual information sources (e.g. spatial edge and geometrical structures) (Moser, Serpico, and Benediktsson 2013). Tarabalka et al. (2010b) integrated the ‘fuzzy no-edge/edge’ function into the energy function to regularize class memberships while preserving the spatial edge information, and demonstrated that MRF-based spatial regulation methods yield a more accurate classification result compared to other classification methods (e.g. SVM, extraction and classification of homogeneous objects, and watershed segmentation).

2.5. Accuracy assessment

Validation, or accuracy assessment of the classification result, is an essential part of land-cover mapping(Congalton et al. 2014).

The accuracy of the classified image has two aspects: geometric accuracy, which is the accuracy related to the geolocation of each pixel or object, and thematic accuracy, which is the accuracy related to the labelling of each pixel or object (Möller et al. 2013; Zhan et al. 2005). Of more interest to users is thematic accuracy, whereas geometric accuracy has a major impact on the representation of thematic accuracy (Foody 2002; Pontius and Millones 2011).

2.5.1. Thematic accuracy

2.5.1.1. Object-based assessment

For pixel-based accuracy assessment, one widely recognized routine is to derive the accuracy metrics from the confusion matrix (CM) (Congalton and Green 2009; Foody 2002). In the case of object-based assessment, the CM can be modified to per-object CM or area-based CM (Lizarazo 2014; Stehman and Wickham 2011).

The most widely used accuracy assessment metric for object-based image classification is overall accuracy, that is, the ratio of the number of correctly classified sampling units to the total number of selected sampling units, which is the per-object accuracy (Radoux et al. 2011; Zhan et al. 2005). However, objects are of varying size and area. For large homogeneous areas (e.g. water), to consider the large object as one unit could result in misleading accuracy estimation. Hence, when the sampling unit is an object, it is common to modify the accuracy metric by using the area of the object, instead of the number of objects, in the calculation, which is the area-weighted or area-based accuracy (Stehman and Czaplewski 1998; Desclée, Bogaert, and Defourny 2006). Another consideration of the object-based accuracy measure is that the object may contain heterogeneity and thus should not be binarized as the ‘correct’ or ‘incorrect’ object (Stehman and Wickham 2011). The theme similarity measure proposed by Lizarazo (2014) considers the ratio of overlapping area that has thematic category matching the area of the validation object.

It should be noted that the accuracy derived from validation samples is not a direct estimate of the accuracy of the entire study area. Therefore, based on the idea of area-weighted accuracy, Radoux et al. (2011) and Radoux and Bogaert (2014) proposed modified predictors of accuracy metrics by considering both the sampled and unsampled objects. The correctness of the unsampled objects was estimated using the statistical expectation. They demonstrated that the proposed predictors are also unbiased estimates of the accuracy metrics, but with improved precision.

2.5.1.2. Multi-scale assessment

There have been efforts to assess the accuracy of image classification results at multiple scales. For example, Wardlow and Callahan (2014) assessed the MODIS irrigated agriculture data set at both state level and field level. State-level assessment provides a general summary of the quality of the classification, while per-field level assessment reveals how the classification error is related to crop type and location. Furthermore, they found that fields of larger size were classified with higher accuracy.

Ge et al. (2009) adopted the rough set theory, which is a useful tool for obtaining knowledge from incomplete or inconsistent data, to assess the thematic accuracy of a classification result in a multi-level fashion without the use of validation data. The three scales considered were pixel-level, object/category-level, and map-level. At the pixel-level, Shannon entropy was used to quantify the amount of information required to label a pixel as a particular class with complete certainty. Rough degree and rough entropy, which measure the uncertainty of labelling the object as a certain class, are used to measure object/category-level accuracy. At the image scale or map level, the accuracy and quality of approximation are used to represent the overall uncertainty. An experiment with five different data sets revealed that object/category-level measures have significantly negative correlation with user’s accuracy and class-specific kappa coefficient, indicating the power of the proposed metrics.

2.5.1.3. Geostatistical error assessment

Thematic error is in many cases not randomly distributed in space, but correlated with the position of boundaries between different classes (Foody 2002). Geostatistical methods are employed to model the spatial pattern of errors. Kyriakidis and Dungan (2001) proposed a method to quantify local classification accuracy by indicator kriging. Another commonly adopted approach is sequential indicator simulation, which investigates the uncertainty propagation in the classification procedure (Magnussen and de Bruin 2003; de Bruin 2000). Carneiro and Pereira (2014) extended geostatistical stochastic simulation to integrate the influence of patch size by considering patch size as one source of uncertainty. On the other hand, geostatistical classification frames are also able to provide uncertainty maps as a by-product (Li and Zhang 2011b).

2.5.2. Geometric accuracy

Geometric accuracy concerns the accuracy associated with space; the labelling of pixels or objects is of no interest at all. For pixel-based classification, geometric accuracy is commonly referred to as positional accuracy (Congalton and Green 2009). For object-based classification, geometric metrics characterizing more properties have been developed, including shape accuracy, edge accuracy, and segmentation accuracy, as well as positional accuracy.

2.5.2.1. Positional accuracy

Pixel-based assessment: due to registration discrepancy, the pixels in the classified image may not coincide perfectly with the validation image or samples. The two-dimensional horizontal accuracy of the classification result can be quantified using the root mean square error (RMSE) between corresponding pixels, which is an overall positional accuracy across the entire image (Congalton and Green 2009). The correspondence of the pixels should be determined with interpreted points from both the classified image and the validation image or samples.

Object-based assessment: area-based metrics quantify how much and how well the classified and validation objects overlap. For example, Persello and Bruzzone (2010) considered the situation where more than one object from the classification map overlaps with one validation object. The classified object that has the largest overlapping area with the validation object is considered to be the corresponding one. Another common approach is location-based, quantifying the deviations of the corresponding object centres. For example, Zhan et al. (2005) proposed a quality measure for the position, Q_Loc, which is the Euclidean distance between the centres of mass of the two corresponding objects. Based on Q_Loc, position similarity indices are designed by normalizing Q_Loc and subtracting this from 1. The normalization term used includes the maximum distance from the furthest validation centroid to the centroid of the overlapping area (Möller et al. 2013), the square root of the overlapping area (8685–8698), the diameter of a combined area circle (CAC, i.e. a circle that has the same area as the sum of the areas of the two corresponding objects) (Lizarazo 2014), and other measures. The RMSE metric can also be applied to object centroids as a global metric to measure the overall positional error of the objects (Whiteside, Maier, and Boggs 2014). The other global metrics include the mean and variance of the accuracy metrics for each object (Zhan et al. 2005).

2.5.2.2. Shape accuracy

Shape accuracy quantifies the similarity between the shape of a validation object and that of the corresponding classified object. The shape of an object can be described using various well established indices, e.g. compactness, eccentricity, or fractal dimension (McGarigal, Cushman, and Ene 2012; Angel, Parent, and Civco 2010). Shape accuracy can be estimated from these shape indices. The dissimilarity of the shape indices for the two corresponding objects is the shape-related error, which is also a measure of the accuracy. Ways of comparing the two shape indices include the absolute value difference (Persello and Bruzzone 2010), the ratio of the values (Lizarazo 2014), and other approaches.

2.5.2.3. Edge accuracy

Edge accuracy quantifies the coincidence between the validation object and the classified object. Usually the ratio of the length of spatially overlapping edges to the length of the validation object edge is calculated to represent the edge accuracy (Persello and Bruzzone 2010; Lizarazo 2014).

Object fate analysis (OFA) was presented by Schöpfer and Lang (2006) and has been applied to evaluate the deviation of object boundaries (Albrecht, Lang, and Hölbling 2010; Hernando et al. 2012). For each validation object, six scenarios are defined according to its spatial relationship with the classified objects. Two metrics, extrution and intrusion, are calculated for each object to quantify the difference in object boundaries (Albrecht, Lang, and Hölbling 2010). The mean value and variance of these metrics can be used to measure the global accuracy (Schöpfer and Lang 2006).

2.5.2.4. Segmentation accuracy

For object-based image classification, segmentation is the procedure for creating the objects. When one validation object is divided into more than one object in the segmentation, over-segmenation occurs. When the segmentation produces objects that encompass more than one validation object, under-segmentation occurs. Over- and under-segmentation are the two major errors associated with segmentation. A large number of metrics for quantifying these two situations have been developed to represent segmentation accuracy. Over-segmentation metrics are mainly based on the ratio of overlapping area to validation object area. Under-segmentation metrics, on the other hand, are mainly based on the ratio of overlapping area to classified object area (Clinton et al. 2010; Persello and Bruzzone 2010). A comprehensive summary of these metrics can be found in Clinton et al. (2010). Zhang et al. (2015a) proposed precision as a global metric for under-segmentation assessment, by taking the ratio of all overlapping areas to all classified object areas; and recall for over-segmentation assessment, by taking the ratio of all overlapping areas to all validation object areas.

Combined metrics have also been proposed to consider both over- and under-segmentation in an integrated fashion (Zhang et al. 2015a; Clinton et al. 2010). Intuitive ways to combine these two metrics include the sum and root mean square (289–299; Möller et al. 2013). The F-measure has also been adopted for this combination (Zhang et al. 2015a).

2.5.3. Combination of geometric and thematic accuracies

Geometric and thematic accuracy measures are typically designed separately and are complementary to each other. Nevertheless, attempts have been made to combine these measures.

Lizarazo (2014) proposed the STEP similarity matrix to comprehensively evaluate the accuracy of object-based image classification. The acronym STEP denotes shape, theme, edge, and position, respectively. The STEP matrix resembles the traditional CM, but shows the four similarities in each cell. Traditional area-based thematic CM, as well as shape CM, edge CM, and position CM, can be derived from the STEP similarity matrix.

Hernando et al. (2012) combined the OFA concept with the traditional pixel-based error matrix concept to create an OFA matrix, to integrally evaluate the thematic accuracy of the classification result in an integrated approach.

Pattern-based accuracy assessment is another approach used to combine geometric and thematic accuracy assessments (Foody 2008). A traditional way of doing pattern-based assessment is to extract pattern indices (e.g. number of target objects, size of target objects) from both classified and validation maps, and to compare the similarity (Taubenbock et al. 2011; Potere et al. 2009). Dihkan et al. (2013) proposed another pattern-based procedure to assess the accuracy of tea plantation mapping. A fuzzy local matching algorithm (Power, Simms, and White 2001) was introduced to measure the local similarity of patterns between the classified and validation maps. A map of the local similarity level is used to show the pattern-based accuracy both quantitatively and qualitatively. Furthermore, global accuracy is calculated from the weighted average of the local similarity measure for each object.

3. Incorporating spatial information in spectral unmixing

Spectral unmixing is the process of decomposing the measured spectrum of a mixed pixel into a set of pure spectral signatures called endmembers and their corresponding abundances, indicating the fractional area coverage of each endmember present in the pixel (Keshava and Mustard 2002). Early research on spectral unmixing dates back three decades (Adams, Smith, and Johnson 1986). Since then, many studies have exploited the spectral information inherent in remotely sensed imagery to analyse mixtures of distinct materials, treating each pixel as independent of its spatial neighbours. Not until recently has spatial information been accounted for in the development of spectral unmixing methods (Van Der Meer 1999; Rand and Keenan 2001; Roessner et al. 2001; Plaza et al. 2002; Maselli 2001). It has been demonstrated that the integration of both spatial and spectral information can lead to improvements in the unmixing results (Shi and Wang 2014). In the remainder of this section, we summarize the state-of-the-art spatial spectral unmixing methods according to the following three categories: 1) endmember extraction, 2) selection of endmember combinations, and 3) abundance estimation.

3.1. Endmember extraction

The technique of endmember extraction aims to determine endmember spectra from the remotely sensed image itself, taking advantage of the fact that image endmembers share the same spatial scale and atmospheric conditions as the image to be unmixed. A substantial number of endmember extraction algorithms have been proposed in the literature (Bioucas-Dias et al. 2012). Based on the theory of convex geometry, many algorithms search for the set of pixels that are located in the vertices of a simplex fit to the convex hull of the pixels in spectral space. However, these geometrically based algorithms utilize spectral properties of the remotely sensed imagery alone and tend to be susceptible to outlier pixels and noise. This problem can generally be alleviated by taking into account the spatial context of the pixels. The incorporation of spatial information in endmember extraction has culminated in the development of spatially oriented algorithms, the refinement of spectral-only algorithms, and the development of preprocessing modules. Spatially oriented algorithms are a fundamentally different approach to geometrically based algorithms – the former rely on detecting pure pixels or pure neighbourhoods in the spatial context, rather than investigating spectral extremeness of pixels in the spectral space. The refined spectral-only algorithms make use of spatial information to improve the performance of geometrically based algorithms; the primary structures of the original spectral-only algorithms are not changed substantially. The preprocessing methods prior to endmember extraction have a unique modularity; no modification of existing endmember extraction algorithms is needed to apply the preprocessing methods. Extracted endmembers that account for the spatial context of candidate pixels may be more spatially representative and less susceptible to outlier pixels.

Spatially oriented endmember extraction algorithms provide a way to assess the purity of a pixel or a spatial neighbourhood without resorting to geometrically based algorithms (Torres-Madronero and Velez-Reyes 2014; Plaza et al. 2002; Mei et al. 2010). The automated morphological endmember extraction (AMEE) algorithm proposed by Plaza et al. (2002) is the earliest attempt in this direction. This algorithm extends conventional greyscale morphology to multi-channel imagery and uses two morphological operations, i.e. dilation and erosion, to identify the respective pixels with highest and lowest purity in a sliding window. A morphological eccentricity index (MEI) value, defined as the spectral angle distance between the purest pixel and the most mixed pixel, is assigned to the purest pixel in the window. Endmembers are obtained by applying a threshold to the MEI image. More recently, Mei et al. (2010) developed the spatial purity-based endmember extraction (SPEE) algorithm to detect pure spatial neighbourhoods. The spatial neighbourhood purity index is determined by either an intensity-level measurement used in the original intensity domain or feature-level measurement used in the transformation domain. A threshold is then applied to obtain endmember candidates. A graph-based spatial refinement process is adopted to reduce the number of endmember candidates. Pure spatial neighbourhood can also be identified using multi-scale representation of remotely sensed imagery (Torres-Madronero and Velez-Reyes 2014). A series of smoothed images are constructed and endmember spectra are extracted from the multi-grid structure. The smoothed spectra are equivalent to the average of spectrally similar and spatially adjacent pixels.

In addition, a collection of endmember extraction algorithms that incorporate spatial information has been proposed to refine spectral-only algorithms (Li and Zhang 2011a; Rogge et al. 2007; Zhang, Rivard, and Rogge 2008; Xu, Du, and Zhang 2014; Xu, Zhang, and Du 2015). These either average spectrally similar and spatially adjacent endmember candidates produced by spectral-only algorithms, or select endmember candidates that are present in homogenous neighbourhoods. The spatial-spectral endmember extraction (SSEE) algorithm (Rogge et al. 2007) is an improvement on pixel purity index (PPI) (Boardman 1993; Boardman, Kruse, and Green 1995). A set of local eigenvectors generated by applying singular value decomposition (SVD) to non-overlapping image subsets is used as projection vectors for PPI. The obtained endmember candidates are updated by a spatial averaging process to average spectrally similar endmember candidates within a sliding window. Zhang, Rivard, and Rogge (2008) improved orthogonal subspace projection (OSP) (Harsanyi and Chang 1994) and vertex component analysis (VCA) (Nascimento and Dias 2005) by averaging the endmember candidates that met both the spectral similarity and spatial adjacency criteria. Li and Zhang (2011a) refined the iterative error analysis (IEA) approach (Neville et al. 1999) by identifying endmember candidates that had sufficient spectrally similar pixels in a local window. Xu, Du, and Zhang (2014) developed the abundance-constrained endmember extraction (ACEE) algorithm and further incorporated spatial information in ACEE. A unique feature of the proposed spatial-spectral information-based ACEE (SSACEE) is that local homogeneity of endmember candidates is determined using estimated abundances instead of pixel spectra. Xu, Zhang, and Du (2015) used PPI as a coarse screening. In each partitioned image subset, endmember candidates identified by PPI were further evaluated according to local homogeneity criteria in order to obtain an endmember set.

In contrast to the incorporation of spatial information in endmember extraction algorithms, a number of studies consider spatial information through a preprocessing step prior to endmember extraction. This preprocessing module can be combined with any endmember extraction algorithm. In the spatial preprocessing (SPP) algorithm (Zortea and Plaza 2009), a spatially derived weighting factor is used to adjust the original spectral signature on a per-pixel basis so that pixels located in spatially homogenous regions will engage a lower number of adjustments than those in spatially heterogeneous regions. As such, when a geometrically based algorithm is applied to the adjusted image, pixels located in spatially homogenous regions are more likely to be identified as endmembers. Many preprocessing methods adopt segmentation or clustering techniques to partition remotely sensed imagery into homogenous regions and use the mean spectra of each region for endmember extraction algorithms (Thompson et al. 2010; Martin and Plaza 2011; Zhang et al. 2014a). The spatial-spectral preprocessing (SSPP) algorithm (Martin and Plaza 2012) provides a more integrated framework to combine both spatial homogeneity and spectral purity at the preprocessing level. Within each region generated by unsupervised clustering, a subset of spatially homogenous and spectrally pure pixels is identified. Spatial and spectral information can also be combined to develop local neighbourhood weights (Liu et al. 2012). A scene is divided into homogenous and transition regions through weight thresholding. Only pixels that fall in homogenous regions are analysed by endmember extraction algorithms. Several preprocessing methods have recently been developed to reduce the size of the original data set while retaining the accuracy of endmember extraction (Rogge et al. 2012; Lopez et al. 2013; Beauchemin 2014). Rogge et al. (2012) applied local endmember extraction to non-overlapping image subsets to constitute a reduced set of pixels. Lopez et al. (2013) selected pixels that were at the spatial edges and spectral extremes of the image. Beauchemin (2014) exploited locally identified pixels that were spectrally similar to the mean spectra of non-overlapping image subsets to define a simplex, and filtered out the pixels located inside this simplex.

A potential drawback of spatial spectral algorithms for endmember extraction is that these algorithms tend to neglect anomalous endmembers (Zortea and Plaza 2009). To overcome this problem, anomaly detection techniques can be integrated with spatial spectral algorithms to extract both homogenous and anomalous endmembers (Mei et al. 2011; Erturk et al. 2014).

3.2. Selection of endmember combinations

Most spectral unmixing studies assume that the set of endmembers is invariant across the entire scene. However, due to endmember variability arising from changing illumination conditions and physical differences in broadly defined endmember classes, it is more realistic to use multiple endmember spectra to represent the same class (Somers et al. 2011; Zare and Ho 2014). Moreover, the actual number of endmembers contained in a single pixel is usually much lower than the number of endmembers present across the entire scene (Roberts et al. 1998). In order to address these two issues, many spectral-only methods have been developed to allow the number and type of endmembers as well as their corresponding spectral signatures to vary on a per-pixel basis (Somers et al. 2011), e.g. the well-known multiple endmember spectral mixture analysis (MESMA) (Roberts et al. 1998). However, these methods usually suffer from spectral confusion and computational complexity induced by numerous endmember combinations. Considering that nearby pixels are likely to share similar endmember combinations, the incorporation of spatial information in the selection of endmember combinations can alleviate spectral confusion and improve computational efficiency. The spatial spectral methods available for selection of endmember combinations can be categorized into two groups: per-pixel methods that allow each pixel to have a specific set of endmembers, and per-field methods that select a distinct endmember combination for each homogenous field partitioned from remotely sensed imagery.

When pure pixels or image endmembers are widely distributed across an image, a per-pixel strategy is usually applied to select for a mixed pixel a set of endmembers that are spatially close to the pixel. Spectrally pure pixels can be identified by supervised classification with rigorous criteria (Roessner et al. 2001; Deng and Wu. 2013; Wu, Deng, and Jia 2014; King and Younan 2006), manual interpretation (Johnson, Tateishi, and Kobayashi 2012; Zhong, Wang, and Wu 2015), maximum abundances of global endmembers (Cui, Li, and Zhao 2014), histogram thresholding (Ma et al. 2014), or purity determination algorithms (Mei, Du, and He 2015). For a mixed pixel located at the centre of a predefined spatial neighbourhood, the optimal per-pixel endmember set can be selected by an iterative unmixing procedure similar to MESMA (Roessner et al. 2001; Wu, Deng, and Jia 2014), an endmember scoring scheme (King and Younan 2006; Mei, Du, and He 2015), or an endmember synthesis process via spatial averaging (Deng and Wu. 2013; Ma et al. 2014; Zhong, Wang, and Wu 2015). It is worth noting that Johnson, Tateishi, and Kobayashi (2012) used a more sophisticated spatial interpolator, i.e. ordinary kriging (Van Der Meer 2012), to synthesize the vegetation indices of bare soil (and/or dead vegetation) and full green vegetation cover, and produced more accurate estimates of fractional green vegetation cover than through the use of inverse distance weighting (IDW) for synthesis. In addition to interpolation, extrapolation has been used to synthesize endmember spectra for coarse-resolution images where pure pixels are absent (Maselli 2001; Zhang et al. 2015b). With known endmember fractions of training samples, multivariate linear regression models are constructed using the fractions as independent variables and the spectra of the sample pixels as dependent variables. Endmember spectra are then synthesized by extrapolating the relevant regression model to 100% cover fraction for each endmember. In order to incorporate endmember variability, Maselli (2001) used locally calibrated regression statistics, including mean vectors and variance/covariance matrices, to construct regression models. Zhang et al. (2015b) adopted geographically weighted regression (GWR) (Brunsdon, Fotheringham, and Charlton 1996) to derive synthetic endmember spectra. Both methods can produce spatially variable endmember signatures. Furthermore, an improvement on MESMA has been made through a hierarchical approach (Franke et al. 2009). Classification results from simpler levels are used as spatial constraints for selection of endmember combinations at more detailed levels. It should be noted that all the aforementioned per-pixel methods represent endmembers as a discrete set. Nonetheless, endmembers can also be represented using a continuous distribution. A recent study exploited the beta compositional model (BCM) for spectral unmixing (Du et al. 2014). Each pixel was considered to follow the distribution of a convex combination of beta-distributed endmembers, which was also approximated by beta distribution. In the proposed BCM-spatial algorithm, K-nearest spatial-spectral neighbours were identified and used to estimate the statistical moments of the beta distribution for each pixel. These statistical moments were then input to the objective function or the likelihood function to estimate the abundances.

The per-field methods for selection of endmember combinations generally partition a heterogeneous scene into a number of homogeneous fields and assign each field a small set of endmembers that are present in the field. A variety of image partitioning approaches have been adopted, including Markov random fields (MRFs) (Rand and Keenan 2001), iterative self-organizing data analysis techniques A (ISODATA) (Shoshany and Svoray 2002), unmixing results (García‐Haro, Sommer, and Kemper 2005), local tiles (Canham et al. 2011; Goenaga et al. 2013), iterative fuzzy clustering (Zare et al. 2013), image classification (Zhang et al. 2014c), and multi-resolution image segmentation (Zhang et al. 2014a; Li et al. 2015). The image partitions so produced are commonly grouped to form homogenous fields that represent different land-cover types according to prior knowledge of the scene (Rand and Keenan 2001; Shoshany and Svoray 2002; García‐Haro, Sommer, and Kemper 2005; Zhang et al. 2014c). This image stratification process is also suggested by a comparative study for estimating fractional green vegetation cover (Xiao and Moody 2005). Endmembers in each field are then selected by a manual process (Rand and Keenan 2001; Shoshany and Svoray 2002; García‐Haro, Sommer, and Kemper 2005; Zhang et al. 2014c; Li et al. 2015) or automatic endmember extraction algorithms (Canham et al. 2011; Goenaga et al. 2013; Zhang et al. 2014a). Zare et al. (2013) proposed a piecewise convex multiple-model endmember detection (PCOMMEND) algorithm to partition an image into multiple fields constrained by a family of simplices in spectral space. A joint optimization procedure was conducted to estimate endmember spectra, fractional abundances, and fuzzy membership values assigned to each pixel denoting the field to which it belonged. An improved algorithm was developed to add spatial constraints on fuzzy membership values in order to guide neighbouring pixels to share the same set of endmembers (Zare et al. 2010).

3.3. Abundance estimation

Fractional abundances of endmembers present in mixed pixels are the primary product of spectral unmixing. The spectrum of a pixel is often modelled as a linear combination of endmember spectra weighted by their corresponding abundances (Settle and Drake 1993; Shimabukuro and Smith 1991), and the abundances can be derived using the least squares method to minimize the mean square error between the actual spectrum and the reconstructed spectrum (Heinz and Chang 2001). In spite of its mathematic simplicity, this method estimates abundances in a pixel-by-pixel fashion and the abundances of each pixel are derived independently of its neighbouring pixels. The incorporation of spatial information, nonetheless, can improve the accuracy of abundance estimation and encourage derived abundances to be more spatially consistent.

An early exploratory work was conducted by Van Der Meer (1999). A series of spatial criteria defined on the RMSE for each pixel were employed for an iterative spectral unmixing (ISU) method with the goal of removing spatial structure remaining in the residuals. However, spatial criteria were not included in an objective function to be optimized. In the spectral and spatial complexity blind source separation (SSCBSS) method (Jia and Qian 2007), the abundance in a pixel was predicted by the weighted average of the abundances in its neighbouring pixels and the objective function was constructed to maximize the spatial predictability of the abundances. A MRF was also exploited for spectral unmixing (Eches, Dobigeon, and Tourneret 2011). In the proposed method, an implicit image classification was conducted to partition the image into homogenous regions where the statistical moments of fractional abundances were unchanged. A hierarchical Bayesian model was then adopted to infer the distributions of the abundances and class labels. In addition, local homogeneity was exploited to select landmark points for nonlinear manifold-based spectral unmixing (Bottai et al. 2013). Recently, image partitioning was integrated with a reweighted λ₁ framework in a spatially weighted sparse coding (SWSC) algorithm (Soltani-Farani and Rabiee 2015).

Arguably, the most popular spatial spectral strategy for abundance estimation makes use of a term of spatial regularization that is involved in the optimization procedure to promote the smoothness of estimated abundances. Examples of regularizers include spatial statistics (Song, Jiang, and Rui 2010), fuzzy local information (Zare 2011), spatial spectral coherence (Castrodad et al. 2011), total variation (Chen, Richard, and Honeine 2014; Iordache, Bioucas-Dias, and Plaza 2012), λ ₂-norm (Liu et al. 2012; Zhu et al. 2014), and others (Jia and Qian 2009; Liu et al. 2011). Spatial regularization is especially important for sparse unmixing, which aims to find the optimal subset of endmember spectra in a spectral library through sparsity-inducing regularization (Iordache, Bioucas-Dias, and Plaza 2011). Since the number of endmembers present in a single pixel is relatively limited, the true abundances associated with a large spectral library are inevitably comprised of many zeros. Therefore, the abundance maps depict pronounced spatial autocorrelation that can be taken into account in the optimization procedure through spatial regularization.

4. Summary and future perspectives

Incorporating spatial information in spectral-based image classification is an important aspect many remote-sensing studies. In this paper, we reviewed and summarized recent advances in various techniques related to this topic and grouped them into five sets: pre-classification methods, sample selection, classifiers, post-classification, and accuracy assessment. Before remote-sensing images are fed to various classifiers, spatial information in most cases is incorporated as additional bands in the original spectral images. GLCM, as a texture measurement, is probably the most widely and frequently adopted method in this case. Recent advancements in other methods to derive texture information, such as fractal dimension, geostatistics, and Gabor and wavelet filters, or image morphological analysis, have also shown their potential to extract spatial information at the pre-classification stage. However, the performance of these different methods may be dependent on the specific application and a comprehensive comparison of these methods has yet to be made.

By including spatial information in sample selection, the sample unit, sample size, and sample design are all affected. A number of studies have incorporated geostatistics into sampling to help determine the sample unit size or even sampling design. However, these studies have mainly focused on interpolation, but not classification. To some extent, to predict the class labels from training samples is also a type of interpolation. Therefore, the conclusions and analytical methods can be borrowed for sample selection in image classification.

GEOBIA, as a new and evolving paradigm in remote sensing and GIScience, provides a spatial perspective for understanding image objects and geographic-based intelligence. In the future, we need to build a more comprehensive geographic knowledge-sharing system with GEOBIA databases, to support collective intelligence and knowledge accumulation. Additionally, object-oriented data models need to be integrated into GEOBIA to more effectively explain the natural and anthropogenic phenomenon with image objects. Two contextual classifiers (SVM and MRF) are primarily reviewed in this study to exploit the spatial information among neighbouring pixels. The introduction of composite kernels in SVM offers valuable opportunities to exploit the spatial information in classifiers. To better take advantage of the diverse spatial information (e.g. shape, size, and homogeneity), more sophisticated spatial kernels need to be explored. MRFs exhibit considerable capacity for incorporating multiple information sources (e.g. multi-scale, multi-bands, multi-temporal) in land-cover mapping. Investigating various information sources in a comprehensive framework is thus needed to advance this classifier.

The emphasis on spatial information has encouraged the development of spatially related accuracy assessment for image classification. Various metrics have been developed for the assessment of both geometric and thematic accuracy of image classification results. Among these, some are conceptually identical and technically similar. A thorough compilation of existing accuracy metrics is needed. In addition, although some interesting topics have been raised, spatial-thematic accuracy assessment is still underdeveloped.

The incorporation of spatial information has also achieved great success for improving the performance of spectral unmixing. In this paper, we summarized the available spatial spectral unmixing methods in the three aspects of spectral unmixing: endmember extraction, selection of endmember combinations, and abundance estimation. Compared with the extent to which spatial information is incorporated in image classification, the area of spatial spectral unmixing is less established and requires more attention from the remote-sensing research community. With this review we hope to trigger new studies to expedite the incorporation of spatial information in spectral unmixing.

Disclosure statement

No potential conflict of interest was reported by the authors.