Spectral-Spatial Diffusion Geometry for Hyperspectral Image Clustering

TL;DR

This paper introduces a spectral-spatial diffusion geometry method for hyperspectral image clustering that leverages spatial regularization to improve clustering accuracy and computational efficiency.

Contribution

It proposes a novel unsupervised clustering algorithm that incorporates spatial regularity into diffusion processes for hyperspectral images, enhancing performance over existing methods.

Findings

Outperforms state-of-the-art clustering algorithms on real hyperspectral data.

Provides a computationally efficient approach with low complexity.

Spatial regularization improves clustering accuracy.

Abstract

An unsupervised learning algorithm to cluster hyperspectral image (HSI) data is proposed that exploits spatially-regularized random walks. Markov diffusions are defined on the space of HSI spectra with transitions constrained to near spatial neighbors. The explicit incorporation of spatial regularity into the diffusion construction leads to smoother random processes that are more adapted for unsupervised machine learning than those based on spectra alone. The regularized diffusion process is subsequently used to embed the high-dimensional HSI into a lower dimensional space through diffusion distances. Cluster modes are computed using density estimation and diffusion distances, and all other points are labeled according to these modes. The proposed method has low computational complexity and performs competitively against state-of-the-art HSI clustering algorithms on real data. In particular, the proposed spatial regularization confers an empirical advantage over non-regularized methods.

Equations

$$d_t(x,y):=\sqrt{{\sum }_{u\in X}(P^t(x,u)-P^t(y,u){)}^2}.$$

$$d_t(x,y)=\sqrt{{\sum }_{n=1}^N{\lambda }_n^{2t}({\Phi }_n(x)-{\Phi }_n(y){)}^2}.$$

$$x\mapsto ({\lambda }_1^t{\Phi }_1(x),{\lambda }_2^t{\Phi }_2(x),\dots ,{\lambda }_M^t{\Phi }_M(x))$$

Full text

I Introduction

As the volume of data captured by remote sensors grows unabated, human capacity for providing labeled training datasets is strained. State-of-the-art supervised learning algorithms (e.g. deep neural networks) require very large training sets to fit the massive number of parameters associated to their high-complexity architectures. In order to take advantage of the deluge of unlabeled remote sensing data available, new methods that are unsupervised—requiring no training data—are necessary.

This article proposes an efficient unsupervised clustering algorithm for hyperspectral imagery (HSI) that exploits not only low-dimensional geometry in the high-dimensional space of spectra, but also the spatial regularity in the 2-dimensional image structure of the pixels. This is achieved by considering spectral-spatial diffusion geometry as captured by spatially regularized diffusion distances [1, 2]. These distances are integrated into the recently proposed diffusion learning algorithm, which has achieved competitive performance versus benchmark and state-of-the-art unsupervised HSI clustering algorithms [3, 4].

The remainder of this article is organized as follows. Background on HSI clustering and diffusion geometry is presented in Sec. II. The proposed algorithm is described and evaluated in Sec. III and IV, respectively. Conclusions and future research directions are presented in Sec. V.

II Background

The unsupervised clustering of HSI $X=\{x_{n}\}_{n=1}^{N}\subset\mathbb{R}^{D}$—understood as a point cloud—consists in providing labels $\{y_{n}\}_{n=1}^{N}$ to each data point without access to labeled training data. The total number of pixels in the image is $N$, and the number of spectral bands is $D$. Typically $D$ is quite large, causing traditional statistical methods to perform sub-optimally. However, the different clusters in the data—corresponding to regions with distinct material properties—often exhibit low-dimensional, though nonlinear, structure. In order to efficiently exploit this structure, methods for clustering HSI that learn the underlying nonlinear geometry have been developed, including methods based on non-negative matrix factorization [5], regularized graph Laplacians [6], angle distances [7], and deep neural networks [8].

Recently, the diffusion geometry [1, 2] of $X$ has been proposed in order to infer the latent clusters. The diffusion distance at time $t$ between $x,y\in X$, denoted $d_{t}(x,y)$, is a notion of distance determined by the underlying geometry of the point cloud $X$. The computation of $d_{t}$ involves constructing a weighted, undirected graph $\mathcal{G}$ with vertices corresponding to the $N$ points in $X$, and weighted edges stored in the $N\times N$ weight matrix $W(x,y):=\exp(-\|x-y\|_{2}^{2}/\sigma^{2})$ if $x\in NN_{k}(y),W(x,y)=0$ otherwise for some scaling parameter $\sigma$ and with $NN_{k}(x)$ the set of $k$-nearest neighbors of $y$ in $X$ with respect to Euclidean distance. Typically $\sigma$ is chose adaptively, for example as some multiple of average distance to nearest neighbors [9], and $k\ll n$ is chosen so that $W$ is sparse; we set $k=100$ in all experiments. Let P(x,y)={W(x,y)}\big{/}{\deg(x)} be a Markov diffusion matrix defined on $X$, where $\deg(x):=\sum_{y\in X}W(x,y)$ is the degree of $x$. The diffusion distance at time $t$ is

[TABLE]

The computation of $d_{t}(x,y)$ involves summing over all paths of length $t$ connecting $x$ to $y$, so $d_{t}(x,y)$ is small if $x,y$ are similar according to $P^{t}$.

The eigendecomposition of $P$ yields fast algorithms to compute $d_{t}$: the matrix $P$ admits a spectral decomposition (under mild conditions, see [2]) with eigenvectors $\{\Phi_{n}\}_{n=1}^{N}$ and eigenvalues $\{\lambda_{n}\}_{n=1}^{N}$, where $1=\lambda_{1}\geq|\lambda_{2}|\geq\dots\geq|\lambda_{N}|$. The diffusion distance (1) can then be written as

[TABLE]

The weighted eigenvectors $\{\lambda_{n}^{t}\Phi_{n}\}_{n=1}^{N}$ are new data-dependent coordinates of $X$, which are nearly geometrically intrinsic [1]. The parameter $t$ measures how long the diffusion process on $\mathcal{G}$ has run: small values of $t$ allow a small amount of diffusion, which may prevent the interesting geometry of $X$ from being discovered. On the other hand, large $t$ allow the diffusion process to run for so long that the fine geometry is homogenized. In all our experiments we set $t=30$; see [4] and [10] for empirical and theoretical analyses of $t$, respectively.

If the underlying graph $\mathcal{G}$ is connected, $|\lambda_{n}|<1$ for $n>1$, so that the sum (2) may approximated by its truncation at some suitable $2\leq M\ll N$. In our experiments, $M$ was set to be the value at which the decay of the eigenvalues $\{\lambda_{n}\}_{n=1}^{N}$ begins to taper [4]. The subset $\{\lambda_{n}^{t}\Phi_{n}\}_{n=1}^{M}$ used in the computation of $d_{t}$ is a dimension-reduced set of diffusion coordinates. Indeed, the mapping

[TABLE]

may be understood as a form of dimension reduction from the ambient space $\mathbb{R}^{D}$ to $\mathbb{R}^{M}$. The truncation also enables us to compute only the first $M=O(1)$ eigenvectors, reducing computational complexity.

Diffusion maps consider the data $X$ as a point cloud in $\mathbb{R}^{D}$. If the data has organizing structure beyond its $D$-dimensional coordinates, this information can be incorporated into the diffusion maps construction by modifying the underlying transition matrix $P$. In the case of HSI, each point is not only a high dimensional spectra, but also a pixel arranged in an image. In particular, the HSI enjoys spatial regularity, in the sense that points in a particular class are likely to have their nearest spatial neighbors being in the same class. The incorporation of spatial information into supervised learning algorithms is known to improve empirical performance for a variety of data sets and methods [11, 12, 13]. In this article, we propose to extend the recently proposed diffusion learning unsupervised clustering framework [4] by directly incorporating spatial information into the underlying diffusion matrix $P$.

III Description of Algorithm

The proposed algorithm first constructs a Markov diffusion matrix, $P$, under the constraint that pixels may only be connected to other pixels that are within some spatial radius $r$. Figure 1 shows how nearest neighbors with and without this spatial constraint differ. The construction of $P$ and the corresponding eigenpairs used to compute the diffusion distances are described in Algorithm 1.

Once the eigenpairs are computed, pairwise diffusion distances are simply Euclidean distances in the new coordinate system (3). Cluster modes are computed as points maximizing $\mathcal{D}_{t}=p(x)\rho_{t}(x)$, where $p(x)$ is a kernel density estimator and $\rho_{t}(x)$ is the diffusion distance of a point to its nearest neighbor of higher density. The mode detection algorithm is summarized in Algorithm 2; see [4] for details.

From the modes, points are labeled iteratively—from highest density to lowest density—according to their nearest spectral neighbor of higher density that has already been labeled, unless it is the case that such a labeling would strongly violate spatial regularity. In that case, points are labeled according to their spatial nearest neighbors. The spectral-spatial labeling scheme is summarized in Algorithm 3, and its crucial parameters and the role of the labeling spatial regularization are discussed at length in [4].

The proposed method—consisting of Algorithms 1, 2, and 3—is called spatially-regularized diffusion learning (SRDL).

IV Experimental Results

We evaluate the SRDL algorithm on $3$ real HSI datasets. The Indian Pines, Salinas A, and Kennedy Space Center datasets considered are standard, have ground truth, and are publicly available111http://www.ehu.eus/ccwintco/index.php?title=Hyperspectral_Remote_Sensing_Scenes. In order to quantitatively evaluate the unsupervised results in the presence of ground truth (GT), we consider three metrics. Overall accuracy (OA) is total number of correctly labeled pixels divided by the total number of pixels, which values large classes more than small classes. Average accuracy (AA) is the average, over classes, of the OA of each class, which values small classes and large classes equally. *Cohen’s *$\kappa$-statistic ($\kappa$) measures agreement across two labelings adjusted for random chance [14].

We note that the Indian Pines and Kennedy Space Center datasets are restricted to subsets in the spatial domain, due to well-documented challenges of unsupervised methods for data containing a large number of classes [15]. These datasets are restricted to reduce the number of classes and achieve meaningful clusters. The Salinas A dataset is considered in its entirety. We remark that results on small subsets can be patched together [4]; such results are not shown here for reasons of space. While the proposed method automatically estimates the number of clusters based on the decay of $\mathcal{D}_{t}$, the number of class labels in the ground truth images were used as parameter $K$ for all clustering algorithms to make a fair comparison with methods that cannot reliably estimate the number of clusters.

Since the proposed and comparison methods are unsupervised, experiments are performed on the entire dataset, including points without ground truth labels. The labels for pixels without ground truth are not accounted for in the quantitative evaluation of the algorithms tested.

IV-A Comparison Methods

We consider 13 benchmark and state-of-the-art methods of HSI clustering for comparison. The benchmark methods are: $K$-means [16] run directly on $X$; principal component analysis (PCA) followed by $K$-means; independent component analysis (ICA) followed by $K$-means[17, 18, 19]222https://www.cs.helsinki.fi/u/ahyvarin/papers/fastica.shtml; Gaussian random projections followed by $K$-means [20]; DBSCAN [21]; spectral clustering (SC) [22]; and Gaussian mixture models (GMM) [23], with parameters determined by expectation maximization.

The recent, state-of-the-art clustering methods considered are: sparse manifold clustering and embedding (SMCE) [24, 25]333http://vision.jhu.edu/code/, which fits the data to low-dimensional, sparse structures, and then applies spectral clustering; hierarchical clustering with non-negative matrix factorization (HNMF) [5]444https://sites.google.com/site/nicolasgillis/code, which has shown excellent performance for HSI clustering when the clusters are generated from a single endmember; a graph-based method based on the Mumford-Shah segmentation [26][6], related to spectral clustering, and called fast Mumford-Shah (FMS) in this article555http://www.ipol.im/pub/art/2017/204/?utm_source=doi; fast search and find of density peaks clustering (FSFDPC) algorithm [27], which has been shown effective in clustering a variety of data sets; and two variants of the recently proposed diffusion learning algorithm, in which the labeling process considers only spectral information (DL) or spectral and spatial information (DLSS) [4].

Among the comparison methods, the proposed method bears closest resemblance to the DL and DLSS methods. SRDL and these two methods differ primarily in how the underlying geometry for clustering is learned. In DL and DLSS, the diffusion geometry is computed though $P$ by considering the HSI only as a spectral point cloud. SRDL regularizes the construction of $P$ by incorporating spatial information into the nearest neighbors construction. The proposed SRDL method also bears similarity to SC, SMCE, and FMS since all these methods use data-driven graphs. The FSFDPC algorithm uses a mode detection scheme similar to SRDL, but without diffusion geometry or spatial information.

IV-B Indian Pines Data

The Indian Pines dataset used for experiments is a subset of the full Indian Pines datasets, consisting of three classes that are difficult to distinguish visually; see Figure 2. Clustering results for Indian Pines are in Figure 3 and Table I. SRDL was run with $r=8$. The spatial regularization in the construction of $P$ is beneficial in terms of labeling: the proposed method improves over DLSS. As expected, SRDL leads to a spatially smoother labeling. However, a mistake is still made in the labeling of the proposed method, indicating that this is a challenging dataset to cluster without supervision.

IV-C Salinas A Data

The Salinas A dataset (see Figure 1) consists of 6 classes in diagonal rows. Certain pixels in the HSI have the same values; in order to distinguish these pixels, small Gaussian noise (variance $<10^{-3})$ was added as a preprocessing step. SRDL was run with $r=20$. Visual results for Salinas A appear in Figure 4 and quantitative results appear in Table I. The proposed method yields the best results, and moreover the labels recovered by the proposed method are quite spatially regular.

IV-D Kennedy Space Center Data

The Kennedy Space Center dataset used for experiments consists of a subset of the original dataset, and contains four classes. Figure 5 shows the data along with ground truth consisting of the examples of four vegetation types which dominate the scene. SRDL was run with $r=20$. Results appear in Table I; for reasons of space, visual results are not shown.

IV-E Parameter Analysis

The crucial parameter in the proposed method is the spatial radius $r$ which determines how near the nearest neighbors in the underlying diffusion process must be. The impact of this parameter in terms of overall accuracy, average accuracy, and $\kappa$ appear in Figure 6. The plots exhibit the tradeoff typical of regularization in machine learning: insufficient or excessive regularization are both detrimental. It is critical to find a good range of regularization parameters, and the flat regions near the maxima in Figure 6 suggest the proposed method is relatively robust to the choice of $r$. We remark that if knowledge of the spatial smoothness of the image is known or can be estimated a priori, then $r$ can be estimated.

V Conclusions and Future Directions

The incorporation of spatial regularity into the construction of the the underlying diffusion geometry improves empirical performance of diffusion learning for HSI clustering in all three datasets considered. For images whose labels are sufficiently smooth, our results suggest there will be a regime for choices of spatial window in which incorporating spatial proximity improves the underlying mode detection and consequent labeling of HSI.

In terms of computational complexity, it suffices to note that the bottleneck is in the construction of $P$. Since $k$ nearest neighbors are sought, and neighbors are constrained to live within a spatial radius $r$, as long as $r,k=O(1)$ with respect to $n$, the nearest neighbor searches for all points can be done in $O(n)$. This gives an overall complexity for the algorithm that is essentially linear in $n$. When the spatial radius is is large enough so that the full scene is considered in the nearest neighbor search, indexing structures (e.g. cover trees) allow for fast nearest neighbor searches, giving a quasilinear algorithm.

The proposed method is likely to be of substantial benefit for scenes that are in some sense smooth with respect to the underlying class labels. For images with many classes that are rapidly varying in space—for example, urban HSI—alternative approaches for incorporating spatial features may be necessary. One approach would be to consider as underlying data points not individual pixels, but higher order semantic features, for example image patches [28]. This would allow for information about fine features such as edges and textures to be incorporated into the diffusion process, allowing for fine-scale clusters to be learned. On the other hand, using patches instead of distinct pixels as the underlying dataset will increase the dimensionality of the data from $D$ to $D\cdot\text{(Patch Size)}$, which may lead to slower computational time. However, it is expected that the patch space is low-dimensional [29], making the statistical learning problem tractable.

Acknowledgments

This research was partially supported by NSF-ATD-1737984, AFOSR FA9550-17-1-0280, and NSF-IIS-1546392.