A Wavelet-Feature Correlation-Based Greedy Clustering Framework for Unsupervised Hyperspectral and Multispectral Remote Sensing Classification

1. Introduction

Hyperspectral remote sensing provides hundreds of contiguous spectral bands, enabling fine discrimination of land-cover materials and supporting applications such as environmental monitoring, precision agriculture, mineral exploration, ecological assessment, disaster monitoring, and military reconnaissance. Supervised deep-learning approaches have recently achieved remarkable performance. Firat et al. developed a hybrid 3D residual spatial–spectral convolutional network that leverages 3D CNN and ResNet architectures to extract deep spatio-spectral features [1,2]. Since ground-truth labels are often unavailable or expensive to acquire, unsupervised clustering remains an indispensable technique for organizing hyperspectral data into spectrally homogeneous regions before further analysis. Unsupervised clustering continues to attract considerable interest because it can exploit large volumes of unlabeled hyperspectral imagery and provide reliable initial information for image interpretation, target detection, and change analysis.

Early hyperspectral image clustering methods mainly relied on distance-based statistical algorithms such as K-means, ISODATA, and fuzzy C-means (FCM). K-means is computationally efficient but requires the number of clusters to be specified in advance [3,4,5] and is sensitive to initialization; its convergence speed is slow [6,7]. The clustering of each pixel reaches a balanced state through iteratively searching for the minimal spatial distance with all class centers. ISODATA introduces automatic cluster splitting and merging but depends heavily on multiple empirical thresholds [8,9,10]. Although ISODATA and ISMC [11,12] can determine the number of classes through self-iteration, the parameters are hard to determine, and some distance parameters should be adjusted according to changes in the number of dimensions. The Orthogonal Subspace Projection (OSP) method [13,14] selects primary targets or directions as a basis. It performs orthogonal projection on all pixels, subsequently identifying suboptimal targets or directions within the projected data. While the primary basis in the OSP method captures the vast majority of information—thereby reducing data dimensionality—the classification results are expressed as grayscale values of mixed pixels relative to the various basis vectors; consequently, the classification performance is not particularly distinct, and the number of spectral bands constrains the number of classes. Cui proposed a feature extraction method by computing vectorized pixel values from a very small local window to yield a superior BoW performance, but this reduces the classification accuracy [15,16]. Peng et al. introduced a graph-based deep spectral–spatial clustering network designed to capture structural relationships among pixels thoroughly. Their approach incorporates a self-expression-embedded multi-graph autoencoder that models high-order structural associations, enabling the extraction of robust spectral–spatial features together with a coherent global clustering structure [17]. FCM improves robustness by assigning fuzzy memberships but often suffers from high computational complexity and sensitivity to noise [18]. Fuzzy C-Means (FCM) improves clustering robustness by assigning soft memberships to pixels; nevertheless, it remains computationally expensive for high-dimensional hyperspectral data and is sensitive to noise and mixed pixels [19]. Graph-based approaches, including normalized cuts and spectral clustering, were subsequently introduced to preserve the intrinsic manifold structure of hyperspectral data and have demonstrated improved clustering performance on complex image scenes [20,21]. More recently, unsupervised representation learning, deep clustering networks, and self-supervised learning have significantly advanced hyperspectral image analysis by learning nonlinear feature representations directly from unlabeled data [22,23,24,25,26].

Although these approaches have achieved encouraging results, several challenges remain. Distance-based clustering methods generally measure similarity using Euclidean distance or related statistical metrics, which often become less discriminative in high-dimensional spectral spaces. Deep clustering approaches usually require computationally intensive network training, large numbers of hyperparameters, and considerable computational resources. Moreover, most existing methods compare complete spectral vectors or learned latent representations while largely overlooking an important characteristic of hyperspectral signatures: although spectral amplitudes may vary substantially because of illumination changes, atmospheric effects, sensor noise, or mixed pixels, spectra belonging to the same material usually exhibit local maxima, local minima, absorption valleys, and abrupt transitions at nearly identical spectral positions. These structural characteristics provide stable descriptors of material properties that are considerably less sensitive to spectral magnitude variations.

Wavelet transforms provide an effective multiscale framework for identifying abrupt spectral variations through zero-crossings and local extrema. Previous studies have applied wavelet analysis to hyperspectral feature extraction and segmentation. Acharyya combined wavelet theory with neuro-fuzzy techniques for segmentation purposes [27,28]. However, most existing wavelet-based methods utilize wavelet coefficients primarily as frequency-domain descriptors without explicitly exploiting the positional consistency of abrupt spectral features for unsupervised clustering. Furthermore, many existing methods still rely on iterative optimization procedures or computationally expensive similarity calculations.

To address these limitations, this paper proposes an Unsupervised Wavelet-Feature Greedy Clustering Algorithm (WFGCA). Instead of measuring similarity directly from spectral amplitudes, the proposed method characterizes each spectrum by multiscale wavelet feature positions corresponding to abrupt spectral structures (Figure 1). A wavelet-feature correlation coefficient is developed to quantify similarity between spectral signatures according to the consistency of these characteristic positions. A greedy clustering strategy is then introduced to sequentially select representative cluster centers with minimal mutual correlation, thereby reducing redundancy while automatically determining the number of clusters without prior knowledge. This greedy clustering provides a global best solution. It resolves the limitations of the previously proposed wavelet-feature correlation ratio Markov clustering algorithm (WFCRMCA) [29], which only provides the best local solutions.

Compared with conventional clustering methods, the proposed WFGCA offers the following contributions: (1) a multiscale wavelet-feature representation that preserves stable spectral structures while reducing sensitivity to spectral magnitude variations; (2) a wavelet-feature correlation coefficient for efficient similarity measurement in high-dimensional hyperspectral data; (3) a greedy clustering algorithm that automatically determines representative cluster centers without requiring the number of clusters to be specified; (4) a computationally efficient unsupervised framework applicable to both multispectral Landsat TM and hyperspectral AVIRIS datasets. Experimental results demonstrate that WFGCA produces accurate clustering results while maintaining low computational complexity and providing flexible control over clustering granularity through several intuitive parameters.

2. Methods

2.1. Wavelet-Feature Correlation Coefficient

A material’s spectral signature is never identical across two acquisitions taken under different conditions. Yet, it tends to rise, fall, peak, or dip at nearly the same wavelength positions each time (Figure 1). This regularity is precisely what the bandpass wavelet kernel ψ(t) [30] is designed to pick up: defined as the first derivative of the Gaussian θ(t) in Equation (1), convolving a spectrum with ψ(t) turns each of the four abrupt-signal shapes in Figure 2a–d, including downward, upward, protruding, and concave, into either a zero-crossing or a local extremum at the filter’s output (Figure 2a′–d′).

$$\theta \left(t\right)={\displaystyle \frac{1}{\sqrt{2\pi }}}{e}^{-t^2/2},\psi \left(t\right)={\displaystyle \frac{d\theta }{dt}}=-{\displaystyle \frac{1}{\sqrt{2\pi }}}t{e}^{-t^2/2}$$

(1)

For some remote sensing data characterized by significant interference or a high proportion of mixed pixels, the locations of abrupt change points may shift, or numerous minor fluctuations may be present. These issues can be addressed by clustering based on selected high-level detail components, applying fuzzy processing to account for positional shifts, or setting a modulus maxima threshold to eliminate insignificant fluctuation signals. Ranked by importance, the four signal types run, from most to least significant: downward, upward, protruding, and concave (Figure 2). In scenes with a large fraction of mixed pixels, the exact location of these features can drift or break into minor spurious fluctuations; setting an upper threshold filters this out, leaving only the strongest minutiae — those that persist at the coarser wavelet scales — for the analysis.

Locating these abrupt signals and recording where they fall is, in essence, what wavelet-feature clustering does. The wavelet coefficients themselves come from Mallat’s fast decomposition algorithm [31,32], Equation (2); from these, upward-step maximum points and downward-step minimum points whose magnitude exceeds a threshold $T_{peak}$ are located and marked.

$$\left\{\begin{array}{c}{c}_{j,k}=\sum _n{h}_{n-2k}{c}_{j-1,n}\\ d_{j,k}=\sum _n{g}_{n-2k}{c}_{j-1,n}\end{array}\right.,j=1,2,\dots ,S_{cale}\left(\le ⌊{log}_{2}b⌋\right)$$

(2)

where $b$ is the total band count, $h$ and $g$ are the low-pass and high-pass filter coefficients, and $c_{j,k}$ and $d_{j,k}$ are the $k$-th low-pass and high-pass filtered values at decomposition level $j$.

Pixel-to-pixel similarity is then quantified by the correlation coefficient (CC), $r_{ij}$, in Equation (3), computed over the $S_{cale2}$-level subset of minutiae taken from the full $S_{cale}$-level decomposition. Rather than Euclidean distance, a geometric measure of how far apart two spectral vectors sit, $r_{ij}$ judges closeness purely from how many partial-minutiae positions the two vectors share, and it is this measure that drives the clustering step.

$$\begin{gathered} r_{ij}={\displaystyle \frac{N\left(t_{i,k}=t_{j,k},k\in \mathsf{\Omega }\right)}{N\left(t_{i,k}\ne t_{j,k},k\in \mathsf{\Omega }\right)+N\left(t_{i,k}=t_{j,k},k\in \mathsf{\Omega }\right)}} \\ if S_{cale2}=S_{cale},\mathsf{\Omega }=\left[0,b{\displaystyle \sum _{k=1}^{S_{cale}}}{\displaystyle \frac{1}{2^k}}-1\right] \\ if S_{cale2}<S_{cale},\mathsf{\Omega }=\left[b{\displaystyle \sum _{k=1}^{S_{cale}-S_{cale2}}}{\displaystyle \frac{1}{2^k}},b{\displaystyle \sum _{k=1}^{S_{cale}}}{\displaystyle \frac{1}{2^k}}-1\right] \end{gathered}$$

(3)

where $t_{i,k}$ is the $k$-th feature value of sampled pixel $i$, and $N(\cdot )$ counts how many feature positions satisfy the stated condition. A binary flag can mark whether a given position is informative enough to be worth including in the comparison. In the special case where $S_{cale2}=S_{cale}=\lfloor {\mathrm{l}\mathrm{o}\mathrm{g}}_{2}b\rfloor$, the number of bits entering the comparison works out to $b{\sum }_{i=1}^{⌊{\mathrm{l}\mathrm{o}\mathrm{g}}_{2}b⌋}1/2^i\approx b-1$.

Working with only a partial set of high-scale minutiae pays off most in scenes with heavy pixel mixing: it keeps the sheer number of minutiae from overwhelming the comparison, at the cost of a mild blurring effect that in practice still produces workable clusters. This same lever, adjusting scale, is what lets the method move between fine-grained and coarse-grained classification. And because comparing abrupt-point locations is computationally cheap, adding more minutiae to the reference set adds little to overall clustering time. As a practical default, expanded multispectral bands are typically processed with $S_{cale2}=S_{cale}$, while hyperspectral bands more often use $S_{cale2}=S_{cale}-2$.

2.2. Expanding Bands Method for Multispectral Images

The number of bands available in multispectral imagery (TM images, for example, provide only seven bands) is often insufficient for the wavelet transform to extract enough informative feature points. To address this, we first exclude bands that are heavily affected by noise, then expand the remaining $b$ bands using second-order and nonlinear correlation functions so that more wavelet features can be detected. This band-expansion method [3] is defined as follows.

• Second-order correlated bands: the auto-correlated band set ${\left\{B_i^2\right\}}_{i=1}^b$ and the cross-correlated band set ${\left\{B_i{B}_j\right\}}_{i,j=1}^b$, $i\ne j$.
• Nonlinearly correlated bands: the square-root band set ${\left\{\sqrt{B_i}\right\}}_{i=1}^b$ and the logarithmic band set ${\left\{\mathrm{l}\mathrm{o}\mathrm{g}{B}_i\right\}}_{i=1}^b$.
• Combining the bands generated in (1) and (2) with the original first-order bands, ${\left\{B_i\right\}}_{i=1}^b$, yields a new dataset containing $(b^2+7b)/2$ bands.

Although this expansion appears to add data-processing complexity, its effect on clustering speed is more favorable than it first seems. When only a few classes are present, each class tends to contain a large number of pixels, and updating the center of a large class is relatively slow. Because the expanded band set resolves more distinguishable classes, individual classes remain smaller on average, which keeps center computation efficient. As a result, the proposed band-expansion method improves both clustering accuracy and clustering speed for multispectral data.

2.3. Greedy Clustering Algorithm

The similarity between two remotely sensed pixels is quantified using the correlation coefficient (CC), which is computed from the number of common abrupt spectral features shared by their spectral signatures. A smaller CC indicates weaker similarity and therefore a larger separation between the two pixels in the feature space. Consequently, the objective of clustering is to identify a set of representative pixels that exhibit the lowest mutual correlations, thereby maximizing the diversity among cluster centers.

The proposed greedy clustering algorithm (GCA) selects representative pixels sequentially according to two criteria. First, a small number of representative samples should capture the dominant information contained in the entire dataset. Second, the correlation coefficient between any pair of selected cluster centers must remain below a predefined threshold, $T_{\mathrm{o}\mathrm{r}\mathrm{t}\mathrm{h}\mathrm{o}}$. This threshold controls the maximum allowable similarity among cluster centers and guarantees that each newly selected center contributes additional information rather than duplicating existing representatives.

Inspired by the concept of orthogonal projection, each selected cluster center is regarded as an independent basis vector describing the feature space. The purpose of selecting a new center is to maximize its independence from the previously chosen centers. Therefore, the cumulative correlation among all selected centers should be minimized. Samples chosen earlier represent the most significant characteristics of the dataset, whereas later selections progressively describe finer details. Because each selection is made based on the current best choice without reconsidering previous decisions, the procedure follows a greedy optimization strategy, referred to as the Greedy Clustering Algorithm (GCA).

The selection process begins with an initial representative sample, denoted as $t_0$. Each subsequent candidate is evaluated according to its correlations with the already selected centers. The optimization objective is expressed by the minimum-energy criterion [33]:

$$\begin{gathered} \mathrm{m}\mathrm{i}\mathrm{n}\mathrm{M}\left(I\right)={\displaystyle \frac{1}{2}}\sum _{j\in I,j\ne i}\sum _{i\in I}{r}_{i,j} \\ \mathrm{s}\mathrm{u}\mathrm{b}\mathrm{j}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{t}\mathrm{o}:r_{i,j}<T_{ortho} \end{gathered}$$

(4)

where $r_{ij}\in \left[\mathrm{0,1}\right]$denotes the correlation coefficient between samples $i$ and $j$; a larger value of $r_{ij}$ indicates stronger similarity between the two samples. The threshold $T_{\mathrm{o}\mathrm{r}\mathrm{t}\mathrm{h}\mathrm{o}}$ plays a role analogous to the error tolerance used in pseudo-inverse computation and determines the admissible correlation among clustering centers.

Following the priority and minimum-energy principles, the GCA evaluates every candidate sample as a potential cluster center. Let $I$ denote the set of selected centers, $J$ the remaining candidate set, and $N_s$ the total number of samples. For each candidate pixel, an importance measure is defined as the sum of its correlation coefficients with all previously selected centers,

$$D_i=\sum _{j\in \mathrm{I}}{r}_{i,j},i\in [0,\dots ,N_s-1]$$

(5)

A smaller value of $D_i$ indicates that the candidate contributes more new information and therefore has a higher priority for selection.

The complete GCA procedure is summarized below.

Initialize the selected set $I=\mathrm{\varnothing }$ and the candidate set $J=\{\mathrm{0,1},\dots ,N_s-1\}$.

Select an initial sample $t_0$ as the first cluster center, move it from $J$ to $I$, and compute the correlation sums $D_i$ for all remaining candidates using (5).

Sort the candidates in ascending order of $D_i$. Beginning with the smallest value, examine each candidate sequentially. A candidate is accepted as a new cluster center only if its correlation coefficient with every previously selected center satisfies

$$r_{ij}<T_{\mathrm{o}\mathrm{r}\mathrm{t}\mathrm{h}\mathrm{o}},\forall i\in I.$$

Once accepted, the sample is transferred from $J$ to $I$, and the correlation sums $D_i$ are updated.

Repeat Step 3 until no additional samples satisfy the correlation constraint. The final set $I$ constitutes the cluster centers.

Assign every pixel to the cluster whose center has the largest correlation coefficient with that pixel.

Unlike conventional iterative clustering approaches, such as K-means, which repeatedly relocate cluster centers until convergence, the proposed GCA directly identifies representative centers through a greedy search. Iterative methods often allocate excessive cluster centers to densely populated regions, increasing redundancy while overlooking subtle but informative structures. Moreover, the required number of clusters must usually be specified in advance.

The proposed GCA avoids these limitations by selecting cluster centers in descending order of information content while enforcing low pairwise correlations. Consequently, the number of clusters is determined automatically according to both the correlation threshold $T_{\mathrm{o}\mathrm{r}\mathrm{t}\mathrm{h}\mathrm{o}}$ and the complexity of the image scene. A larger value of $T_{\mathrm{o}\mathrm{r}\mathrm{t}\mathrm{h}\mathrm{o}}$ imposes a stricter independence requirement, allowing more detailed image structures to be represented and generally producing a larger number of cluster centers. Conversely, scenes with relatively homogeneous land-cover types exhibit higher spectral similarity among pixels, resulting in fewer selected centers.

2.4. Wavelet-Feature Greedy Clustering Algorithm (WFGCA)

The proposed Greedy Clustering Algorithm (GCA) can be extended from the original spectral space to the wavelet feature space by utilizing abrupt spectral features extracted through multiscale wavelet analysis. This extension is referred to as the Wavelet-Feature Greedy Clustering Algorithm (WFGCA). Instead of directly comparing spectral vectors, WFGCA measures pixel similarity based on the distribution of their wavelet-derived feature points, enabling more reliable discrimination of hyperspectral signatures. The implementation of WFGCA is summarized as follows.

• Input parameters: $x_{step}$ and $y_{step}$ are the sampling steps along the horizontal and vertical directions in the remotely sensed images. b is the spectral band number, m is the column number, and n is the row number. $S_{cale}$ is the maximal scale of the applied wavelet transform, normally no more than $⌊{log}_{2}b⌋$. The wavelet minutiae from $S_{cale}-S_{cale2}$ to $S_{cale}$ will be clustered according to (3).
• Apply bandpass wavelet filters [34,35] in (6) to all sampled data at $S_{cale}$ scales, find the positions of abrupt points at each minutiae scale, and mark respectively the upward extreme point and the downward extreme point, whose absolute value is larger than $T_{peak}$, as 1 and 2.

  $$\begin{gathered} \mathrm{H}=\left\{\mathrm{0.0,0.125,0.375,0.375,0.125,0.000}\right\} \\ \mathrm{G}=\left\{-0.0061,-0.0869,-\mathrm{0.5798,0.5798,0.0869,0.0061}\right\} \end{gathered}$$

  (6)
• Sample all the pixels evenly at spatial steps with $x_{step}$ as the x-coordinate step and $x_{step}$ as the y-coordinate step, creating N_s sampled pixels.
• Apply GCA to the minutiae sections for the $S_{cale}-S_{cale2}~S_{cale}$ wavelet coefficients and achieve clustering centers with CC between each other less than $T_{ortho}$.
• Select sample t₀ as the first object, according to (7), initialize the selected sample set I, and the waiting sample set J.

  $$\begin{gathered} I=\left\{t_0\right\},J=\left\{\mathrm{0,1},\dots ,N_s-1\right\}-\left\{t_0\right\} \\ D=\left\{r_{t_{0}j}|j\in \left[0,N_s-1\right]\right\},t_0\in \left[0,N_s-1\right] \end{gathered}$$

  (7)
• In J, choose the sample $i\in \mathit{J}$, with which the CC of each sample in I is less than$T_{ortho}$, and the sum of these CCs is minimal according to (8). If none is found, then clustering ends and goes to step 7; otherwise, add i to object set I and iterate this process.

  $$\begin{gathered} \exists i\in \mathit{J},let{D}_i<D_l\left(i.e.\sum _{j\in I}{r}_{ij}<\sum _{j\in I}{r}_{lj}\right),j\in \mathit{I},l\in \mathit{J}, and r_{lj}<T_{ortho},r_{ij}<T_{ortho} \\ \Rightarrow \mathit{I}=\mathit{I}\cup \left\{i\right\},\mathit{J}=\mathit{J}/\left\{i\right\},D_j=D_j+r_{ij},j\in \left[0,N_s-1\right] \end{gathered}$$

  (8)
• Each pixel $s_i,i\in [0,m*n-1]$, according to the clustering centers created in step 4, is clustered into the class whose center has the maximal CC with $s_i$.

  $$\forall i\in \left[0,m\times n-1\right],\exists k\in \mathbf{I},r_{ij}<r_{ik},j\ne k,j\in \mathbf{I}\Rightarrow s_i=k$$

  (9)

The proposed WFGCA possesses several computational and practical advantages.

• Linear computational complexity. Assuming that $k$ cluster centers are generated from $n$ sampled pixels, only $k(n-1)$ correlation evaluations are required. Consequently, the overall computational complexity is approximately $O\left(n\right)$, making the algorithm suitable for large hyperspectral datasets.
• Efficient similarity computation. Unlike Euclidean-distance-based methods that require extensive floating-point operations, the proposed correlation coefficient is obtained simply by counting matched abrupt spectral features occurring at identical spectral locations. This significantly reduces computational overhead and accelerates the clustering process.
• Low memory consumption. During the greedy search, only the correlation values between candidate pixels and the currently selected cluster centers are maintained. Therefore, the algorithm avoids storing a complete pairwise similarity matrix, resulting in substantially lower memory requirements.
• Adaptive clustering resolution. The clustering granularity can be adjusted by modifying three parameters: $T_{\mathrm{o}\mathrm{r}\mathrm{t}\mathrm{h}\mathrm{o}}$, $S_{cale}$, and $S_{cale2}$. Increasing $T_{\mathrm{o}\mathrm{r}\mathrm{t}\mathrm{h}\mathrm{o}}$ permits finer distinctions among similar spectral patterns, thereby producing a larger number of cluster centers. Conversely, lower threshold values generate fewer clusters by emphasizing only the dominant spectral structures.

3. Results

The characteristics of wavelet-feature clustering were tested using WFMCA analysis of TM and AVIRIS data (Figure 3a and Figure 9a), with classified pixels shown in white.

3.1. AVIRIS Data

Evaluation data were collected by the AVIRIS airborne imaging spectrometer on 2 September 1993, over the SEIDAM Sooke Lake study area on Vancouver Island, British Columbia, Canada [36,37] (Figure 3a). The sensor recorded 224 contiguous spectral bands covering the 410–2450 nm range. Although the sensor’s native ground sampling distance was approximately 20 m, the version of the dataset used here was supplied as a georeferenced product resampled to 25 m. The extracted area has a width of 256 pixels and height of 256 pixels. The WFGCA uses all 224 bands for classification algorithm analysis.

After 5*5 down-sampling, the top three-scale wavelet features from the total five-scale wavelet decomposition ($S_{cale2}/S_{cale}$=3/5) were selected for clustering. When $T_{ortho}$=0.8, 125 classes were created, and the major five clustering results were shown in Figure 3b-5f. The division results are quite ideal. From Table 1, when the 1st and 2nd scale minutiae are used for clustering, the number of classes increases sharply, but clustering accuracy improves only slightly because $T_{ortho}$ controls it. If only the minutiae from the 3rd to 5th scales are used for clustering, and $T_{ortho}$ is increased to 0.95, the number of classes doubles, while clustering accuracy and speed are clearly improved. Therefore, for hyperspectral clustering, partial high-scale minutiae sections and higher T_ortho values can ensure not only a reasonable number of classes but also fairly high clustering accuracy. All the pixels can be identified, according to their spectra, to be forest (Figure 3b), levees (Figure 3c), wetlands (Figure 3d), mountain peaks (Figure 3e), or water (Figure 3f).

3.1.1. Sampling Interval Effect

Increasing the spatial sampling interval x_step sharply reduces both the number of clusters and the WFGCA runtime. Figure 4 and Table 2 show that from x_step = 1 to 5, the number of classes falls from 409 to 125, while runtime decreases from 2652 to 6 s. At x_step = 10, only 81 classes are generated in 3 s, and at 20 the result is 44 classes in 2 s. Thus, full-resolution sampling provides the finest partition but incurs a very high computational cost.

Figure 3. (a) The 50th band Sook Lake AVIRIS image after gray balance. (b)-(f) are the first five major signals after WFGCA clustering with x_step=y_step=5, T_peak=5, T_ortho=0.8, $S_{cale2}/S_{cale}$ = 3/5.

Figure 3. (a) The 50th band Sook Lake AVIRIS image after gray balance. (b)-(f) are the first five major signals after WFGCA clustering with x_step=y_step=5, T_peak=5, T_ortho=0.8, $S_{cale2}/S_{cale}$ = 3/5.

The decline is strongly nonlinear. Increasing x_step from 1 to 2 alone reduces runtime by 95.8% (2652 to 112 s), whereas the class number decreases by 33.5% (409 to 272). Beyond approximately x_step = 9, runtime remains between 2 and 4 s, and the class number fluctuates within a narrower range of 44–81, indicating diminishing computational returns from further subsampling. These small class-count fluctuations reflect changes in which the sampling grid retains representative pixels.

These results reveal a tradeoff between clustering granularity and efficiency: smaller intervals retain more candidate spectra and produce more classes, whereas larger intervals provide faster but coarser partitions. For this AVIRIS dataset, x_step = 5 offers a practical compromise, producing 125 classes in 6 s and reducing runtime by approximately 99.8% relative to x_step = 1.

Figure 4. Effect of $x_{step}$ on (a) class number and (b) time is investigated using WFGCA clustering with x_step=y_step, T_peak=5, T_ortho=0.8, and $S_{cale2}/S_{cale}$ = 5/7.

Figure 4. Effect of $x_{step}$ on (a) class number and (b) time is investigated using WFGCA clustering with x_step=y_step, T_peak=5, T_ortho=0.8, and $S_{cale2}/S_{cale}$ = 5/7.

3.1.2. Band Number

Band number b controls the spectral dimensionality available to WFGCA. Figure 5 shows that both class number and runtime increase nonlinearly as the band number grows. Across the tested range, the class number rises from approximately 1 at b = 2 to nearly 5000 at b = 175, while runtime increases from approximately 1 s to about 1500 s. Growth is relatively gradual at lower band counts but accelerates markedly beyond roughly 125–135 bands, with local fluctuations reflecting differences in the information contributed by individual bands. Additional bands therefore enable finer spectral separation but also expand the feature space and computational burden; band selection or dimensionality reduction may be beneficial when the full spectral resolution is unnecessary.

Figure 5. The effect of band number over (a) class number and (b) time is investigated using WFGCA clustering with x_step=y_step=5, T_peak=5, T_ortho=0.8, and $S_{cale2}/S_{cale}$ = 5/7.

Figure 5. The effect of band number over (a) class number and (b) time is investigated using WFGCA clustering with x_step=y_step=5, T_peak=5, T_ortho=0.8, and $S_{cale2}/S_{cale}$ = 5/7.

3.1.3. T_peak Effect

Increasing T_peak produces a nonlinear reduction in both class number and runtime. In Figure 6 and Table 2, from T_peak = 1 to 5, the class number decreases from 345 to 125, and runtime falls from 14 to 6 s. The decline then slows: at T_peak = 8, WFGCA generates 93 classes in 4 s, and over T_peak = 8–16 the result remains within 83–93 classes and 4–5 s. This plateau indicates that thresholds above approximately 8 provide little additional computational benefit. A low T_peak retains weaker peaks and yields a finer partition, whereas a high value suppresses them and produces fewer classes. For this dataset, T_peak = 5 provides a practical balance, retaining 125 classes with a runtime of 6 s.

3.1.4. T_ortho Effect

T_ortho has a nonlinear effect on clustering granularity. In Figure 7 and Table 3, the class number rises only from 3 to 13 as T_ortho increases from 0.1 to 0.5, then sharply to 28, 59, 125, and 270 at 0.6–0.9, before leveling off at 280 for 0.95. Runtime remains 1–2 s through 0.6 and then rises to 4, 6, 11, and 10 s at 0.7–0.95. Despite small timing fluctuations, computational costs generally increase with class count. Thus, higher T_ortho values produce finer clustering at greater cost, while 0.7–0.8 offers a practical balance for this dataset.

3.1.5. S_cale2 Effect

S_cale2 has a strongly nonlinear effect on clustering granularity and computational cost. From Figure 8 and Table 4, for S_cale2 = 1–3, WFGCA produces only 1, 1, and 2 classes, respectively, with a runtime of 2 s. The class number then increases to 25 at S_cale2 = 4, 125 at 5, 414 at 6, and 2453 at 7, while runtime rises from 1 s at 4 to 6, 26, and 305 s at 5–7. Relative to S_cale2 = 5, including the second- and first-scale minutiae at values of 6 and 7 increases the class number by approximately 3.3 and 19.6 times, respectively, but increases runtime by approximately 4.3 and 50.8 times. Thus, intermediate values of 4–5 provide a more controlled balance, whereas 6–7 should be used only when very fine class separation justifies the substantially greater computational cost.

3.2. TM Data

To evaluate the effectiveness of the proposed Wavelet-Feature Greedy Clustering Algorithm (WFGCA) on multispectral remote sensing imagery, experiments were conducted using the Landsat Thematic Mapper (TM) dataset for the Mississippi River (scene Y4 018 116 055X0 [38], 512*512, 8 bit/pixel) shown in Figure 9a. The test scene contains several representative land-cover types, including water bodies, agricultural fields, vegetation, bare soil, and urban areas. Compared with hyperspectral imagery, Landsat TM provides only a limited number of spectral bands; therefore, discrimination of spectrally similar materials depends primarily on the clustering algorithm’s ability to exploit the characteristic spectral structures contained in the available bands. In WFGCA, the seven bands were expanded to 49 bands, including original data (1-7), second-order auto-correlated bands (8-14), second-order cross-correlated bands (15-35), square root function bands (36-42), and logarithmic function bands (43-49).

The quantitative evaluation summarized in Table 5 further confirms the superiority of the proposed algorithm. After 5*5 down-sampling, the increase in class number from second-order auto-correlation (8-14) is limited to 2, but the expanded bands from second-order cross-correlation (15-35) provide 57 classes and the square root function 90 classes, enhancing classification accuracy. The logarithmic function does not contribute to more detailed classes, even lower down the class number from 90 to 87. Figure 9b-f are the five major classes chosen from the total 87 classes, where the band number is 49, T_ortho is 5, T_ortho is 0.8, $S_{cale}$ uses 5, and $S_{cale2}$ is the same as $S_{cale}$. All the pixels can be identified by their spectra as beach (Figure 9b), plow land or meadow (Figure 9c), river channel (Figure 9d), dyke (Figure 9e), and the slope on the bank (Figure 9f).

Table 5. Effect of $S_{cale2}$ with x_step=y_step=5, T_peak=5, T_ortho=0.8, and $S_{cale}$=7.

${\mathit{S}}_{\mathit{c}\mathit{a}\mathit{l}\mathit{e}2}$		Class number	Time
1		1	2
2		1	2
3		2	2
4		25	1
5		125	6
6		414	26
7		2453	305

Table 5. Effect of $S_{cale2}$ with x_step=y_step=5, T_peak=5, T_ortho=0.8, and $S_{cale}$=7.

${\mathit{S}}_{\mathit{c}\mathit{a}\mathit{l}\mathit{e}2}$		Class number	Time
1		1	2
2		1	2
3		2	2
4		25	1
5		125	6
6		414	26
7		2453	305

Table 6. Mississippi TM image: influence of expanded band number on WFGCA (5*5 sampling).

Band no.	Tpeak	Tortho	Scale	Scale2	Class number	time
7	5	0.8	2	2	1	1
14	5	0.8	3	3	3	14
35	5	0.8	5	5	57	13
42	5	0.8	5	5	90	20
49	5	0.8	5	5	87	20

Table 6. Mississippi TM image: influence of expanded band number on WFGCA (5*5 sampling).

Band no.	Tpeak	Tortho	Scale	Scale2	Class number	time
7	5	0.8	2	2	1	1
14	5	0.8	3	3	3	14
35	5	0.8	5	5	57	13
42	5	0.8	5	5	90	20
49	5	0.8	5	5	87	20

Figure 9. Mississippi TM remotely sensed image: Figure 1. (a) The 3rd band Mississippi TM image after gray balance. b-f are WFGCA clustering the first five main signals with T_ortho=0.8.

Figure 9. Mississippi TM remotely sensed image: Figure 1. (a) The 3rd band Mississippi TM image after gray balance. b-f are WFGCA clustering the first five main signals with T_ortho=0.8.

5. Discussion

The experimental results support the central premise of WFGCA: the positional consistency of multiscale spectral features can provide a compact basis for unsupervised classification when absolute spectral amplitudes vary. In the AVIRIS experiment, the method separated spectra associated with forest, levees, wetlands, mountain peaks, and water. In contrast, the Landsat TM experiment produced spatially coherent groups associated with beach, plowed land or meadow, river channel, dike, and bank slope. These outcomes are consistent with the motivation presented in the Introduction and with the method described in Section 2, because cluster assignment depends on shared extrema and abrupt-transition positions rather than on the Euclidean separation of complete spectral vectors. The findings therefore indicate that wavelet-feature locations retain material-discriminative information in both hyperspectral and multispectral settings.

The AVIRIS parameter analyses clarify how each component of WFGCA controls the balance between clustering granularity and computational cost—increasing the spatial sampling interval from 1 to 5 reduced runtime from 2652 to 6 s (approximately 99.8%) while decreasing the number of clusters from 409 to 125, demonstrating that representative spatial sampling is the main source of computational savings. Increasing T_peak removed weak wavelet extrema and reduced both class count and runtime, with limited additional benefit above approximately 8. T_ortho controlled the admissible similarity among selected centers: values of 0.7-0.8 provided a practical intermediate resolution, whereas larger values produced many more centers at greater cost. S_cale2 had the strongest nonlinear effect; incorporating finer-scale minutiae increased the number of classes from 125 at S_cale2 = 5 to 414 and 2453 at S_cale2 = 6 and 7, respectively, while runtime increased from 6 to 26 and 305 s. Collectively, these results confirm that sampling interval, feature threshold, orthogonality threshold, and decomposition scale serve distinct but complementary roles in controlling efficiency and resolution.

The multispectral experiment further demonstrates why band expansion is needed when only a small number of original bands are available. Expanding the seven Landsat TM bands to 49 bands created additional wavelet-feature positions that the correlation measure could use. However, the contributions of the transformations were unequal: second-order autocorrelation added only two clusters, whereas cross-correlation increased the class count to 57 and the square-root transformation increased it to 90. Adding logarithmic bands did not improve granularity and reduced the final count from 90 to 87. This result suggests that band expansion should be selective rather than automatic; transformations that introduce complementary spectral structure are useful, whereas redundant transformations may add computation without improving discrimination.

Relative to K-means, ISODATA/ISMC, and projection-based approaches discussed in the Introduction, WFGCA offers two practical advantages. First, it does not require a predefined number of clusters; the number emerges from scene complexity and the center-correlation threshold. Second, it avoids repeated center relocation and uses compact feature matching rather than full-vector Euclidean-distance calculations. The resulting clusters are therefore interpretable as groups organized around representative spectra with limited mutual correlation. Nevertheless, the present experiments do not establish statistical superiority over these methods, because no common ground-truth labels or quantitative external validity measures were used. The observed homogeneous regions, preserved boundaries, class counts, and runtimes demonstrate feasibility and parameter behavior, but claims about classification accuracy should be confirmed through controlled benchmark comparisons.

Several limitations should guide interpretation of the results. WFGCA is unsupervised, so its clusters are spectral groupings rather than guaranteed semantic land-cover classes, and the class count can substantially exceed the number of broad cover categories in a scene. The selected parameters remain dataset dependent, particularly at fine wavelet scales where noise and mixed pixels can generate many additional features. Uniform spatial sampling may also omit rare materials that do not fall on the sampling grid, and the current spectral-only formulation does not explicitly use neighborhood context, which could improve robustness in heterogeneous regions.

Future work should therefore evaluate WFGCA on labeled hyperspectral benchmarks using clustering accuracy, normalized mutual information, adjusted Rand index, boundary preservation, runtime, and memory consumption under a common experimental protocol. Adaptive selection of T_peak, T_ortho, sampling interval, and S_cale2 could reduce manual tuning, while rare-target-aware or multistage sampling could limit the loss of small classes. Incorporating spatial-context constraints, automated selection of informative expanded bands, and parallel implementation would further improve robustness and scalability. Extensions to multisensor fusion and temporal remote sensing should be investigated only after these benchmark evaluations establish when positional wavelet features provide advantages over amplitude-based and learned representations.

6. Conclusions

The proposed Wavelet-Feature Greedy Clustering Algorithm (WFGCA) for unsupervised multispectral and hyperspectral remote sensing image classification addresses several limitations of conventional clustering algorithms, including their dependence on distance-based similarity measures, empirical parameter selection, predefined numbers of clusters, and sensitivity to spectral amplitude variations. Instead of comparing complete spectral vectors, WFGCA characterizes spectral signatures using multiscale wavelet-derived feature points and evaluates spectral similarity through a wavelet-feature correlation coefficient that emphasizes the positional consistency of characteristic spectral structures. Representative cluster centers are then selected using a greedy strategy that automatically determines the clustering structure while reducing redundancy among representative spectra.

Experimental evaluations using Landsat TM multispectral imagery and AVIRIS hyperspectral imagery demonstrate the effectiveness of the proposed approach. The parameter analyses show that the spatial sampling interval, wavelet-feature threshold, orthogonality threshold, and wavelet decomposition scale provide effective control over clustering resolution and computational efficiency. Appropriate parameter selection substantially reduces computational cost while preserving the principal spectral characteristics required for accurate clustering. These results indicate that representing spectra by stable wavelet features improves robustness to spectral variability while maintaining the ability to discriminate different surface materials.

An important advantage of WFGCA is its computational efficiency. The combination of spatial sampling, compact wavelet-feature representations, and greedy representative-spectrum selection significantly reduces the number of similarity calculations required during clustering. Unlike iterative optimization algorithms that require repeated updates of cluster centers or empirical splitting and merging operations, the proposed framework automatically identifies representative cluster centers without requiring a predefined number of classes, making it well suited for large-scale remote sensing image analysis.

Although the proposed method demonstrates strong performance on both multispectral and hyperspectral datasets, several opportunities remain for further improvement. Future work will investigate adaptive parameter optimization to reduce manual parameter tuning, incorporate spatial-context information to improve the discrimination of complex heterogeneous regions, and develop parallel implementations suitable for high-performance computing platforms. In addition, extending the wavelet-feature correlation framework to emerging remote sensing modalities, such as multisensor data fusion and time-series image analysis, may further improve the applicability of the proposed algorithm to large-scale Earth observation tasks.

Overall, the proposed WFGCA provides an efficient, robust, and adaptive framework for unsupervised remote sensing image classification. By combining multiscale wavelet-feature extraction, wavelet-feature correlation analysis, and greedy representative-spectrum selection, the method effectively balances clustering accuracy, computational efficiency, and automation, offering a practical solution for the unsupervised analysis of multispectral and hyperspectral remote sensing imagery.