Spectral Algorithms

SPy implements various algorithms for dimensionality reduction and supervised & unsupervised classification. Some of these algorithms are computationally burdensome and require iterative access to image data. These algorithms will almost always execute significantly faster if the image data is loaded into memory. Therefore, if you use these algorithms and have sufficient computer memory, you should load your image into memory.

In [1]: from spectral import *

In [2]: img = open_image('92AV3C.lan').load()

Unsupervised Classification

Unsupervised classification algorithms divide image pixels into groups based on spectral similarity of the pixels without using any prior knowledge of the spectral classes. SPy provides several clustering algorithms for unsupervised classification.

Single-Pass Clustering

The cluster function performs a relatively fast single-pass clustering algorithm. Each pixel in the image is accessed only once (hence “single-pass”), which can make this function useful for a very large image that will not fit into system memory. You must select the number of clusters to generate. cluster will return a classification map and an array containing the locations of the clusters.

In [3]: (m, c) = cluster(img, 20)
Clustering image...done

In [4]: w = view_indexed(m)

Single-Pass clustering results

While cluster generates a set of clusters with a single pass through the data, the clusters are not very accurate and results can be highly dependent on the initial cluster centers used by the algorithm.

k-means Clustering

The k-means algorithm takes an iterative approach to generating clusters. The parameter k specifies the desired number of clusters to generate. The algorithm begins with an initial set of cluster centers (e.g., results from cluster). Each pixel in the image is then assigned to the nearest cluster center (using distance in N-space as the distance metric) and each cluster center is then recomputed as the centroid of all pixels assigned to the cluster. This process repeats until a desired stopping criterion is reached (e.g., max number of iterations). To run the k-means algorithm on the image and create 20 clusters, using a maximum of 50 iterations, call kmeans as follows.

In [5]: (m, c) = kmeans(img, 20, 30)
Initializing clusters along diagonal of N-dimensional bounding box.
Starting iterations.
Iteration 1...done
        21024 pixels reassigned.
Iteration 2...done
        11214 pixels reassigned.
Iteration 3...done
        4726 pixels reassigned.
---// snip //---
Iteration 28...done
        248 pixels reassigned.
Iteration 29...done
        215 pixels reassigned.
Iteration 30...done
        241 pixels reassigned.
^CIteration 31... 15.9%KeyboardInterrupt: Returning clusters from previous iteration

Note that I interrupted the algorithm with a keyboard interrupt (CTRL-C) after the 30th iteration since there were only about a hundred pixels migrating between clusters at that point. kmeans catches the KeyboardInterrupt exception and returns the clusters generated at the end of the previous iteration. If you are running the algorithm interactively, this feature allows you to set the max number of iterations to an arbitrarily high number and then stop the algorithm when the clusters have converged to an acceptable level. If you happen to set the max number of iterations too small (many pixels are still migrating at the end of the final iteration), you can simply call kmeans again to resume processing by passing the cluster centers generated by the previous call as the optional start_clusters argument to the function.

k-means clustering results

Notice that kmeans appears to have captured much more of the spectral variation in the image than the single-pass cluster function. Let’s also take a look at the the cluster centers produced by the algorithm:

In [6]: import pylab

In [7]: pylab.figure()

In [8]: pylab.hold(1)

In [9]: for i in range(c.shape[0]):
   ...:      pylab.plot(c[i])
   ...: 

In [10]: pylab.show()

k-means cluster centers

Supervised Classification

Training Data

Performing supervised classification requires training a classifier with training data that associates samples with particular training classes. To assign class labels to pixels in an image having M rows and N columns, you must provide an MxN integer-valued ground truth array whose elements are indices for the corresponding training classes. A value of 0 in the ground truth array indicate an unlabeled pixel (the pixel is not associated with a training class).

The following commands will load and display the ground truth map for our sample image:

In [11]: gt = open_image('92AV3GT.GIS').read_band(0)

In [12]: w = view_indexed(gt)

Image ground truth classes

We can now create a TrainingClassSet object by calling create_training_classes:

In [13]: classes = create_training_classes(img, gt)

With the training classes defined, we can then create a supervised classifier, train it with the training classes, and then classify an image.

Gaussian Maximum Likelihood Classification

In this case, we’ll perform Gaussian Maximum Likelihood Classification (GMLC), so let’s create the appropriate classifier.

In [14]: gmlc = GaussianClassifier(classes)
  Omitting class   1 : only 54 samples present
  Omitting class   7 : only 26 samples present
  Omitting class   9 : only 20 samples present
  Omitting class  13 : only 212 samples present
  Omitting class  16 : only 95 samples present
Covariance.....done
Covariance.....done
Covariance.....done
Covariance.....done
Covariance.....done
Covariance.....done
Covariance.....done
Covariance.....done
Covariance.....done
Covariance.....done
Covariance.....done

When we created the classifier, it was automatically trained on the training sets we provided. Notice that the classifier ignored five of the training classes. GMLC requires computing the inverse of the covariance matrix for each training class. Since our sample image contains 220 spectral bands, classes with fewer than 220 samples will have singular covariance matrices, for which we can’t compute the inverse.

Once the classifier is trained, we can use it to classify an image having the same spectral bands as the training set. Let’s classify our training image and display the resulting classification map.

In [15]: clmap = gmlc.classify_image(img)
Classifying image...done

In [16]: w = view_indexed(clmap)

Gaussian Maximum Likelihood classification results

The classification map above shows classification results for the entire image. To view results for only the ground truth pixels we must mask out all the pixels not associated with a training class.

In [17]: gtresults = clmap * (gt != 0)

In [18]: w = view_indexed(gtresults)

GMLC results for training set

If the classification results are good, we expect the classification map above to look very similar to the original ground truth map. To view only the errors, we must mask out all elements in gtResults that do not match the ground truth image.

In [19]: gterrors = gtresults * (gtresults != gt)

In [20]: w = view_indexed(gterrors)

GMLC errors

The five contiguous regions in the error image above correspond to the ground truth classes that the GaussianClassifier ignored because they had too few samples.

The following table lists the supervised classifiers SPy currently provides.

Classifier	Description
GaussianClassifier	Gaussian Maximum Likelihood
MahalanobisDistanceClassifier	Mahalanobis Distance
PerceptronClassifier	Multi-Layer Perceptron

Dimensionality Reduction

Processing hyperspectral images with hundreds of bands can be computationally burdensome and classification accuracy may suffer due to the so-called “curse of dimensionality”. To mitigate these problems, it is often desirable to reduce the dimensionality of the data.

Principal Components

Many of the bands within hyperspectral images are often strongly correlated. The principal components transformation represents a linear transformation of the original image bands to a set of new, uncorrelated features. These new features correspond to the eigenvectors of the image covariance matrix, where the associated eigenvalue represents the variance in the direction of the eigenvector. A very large percentage of the image variance can be captured in a relatively small number of principal components (compared to the original number of bands) .

The SPy function PrincipalComponents. This object also contains a transform to rotate data in to the space of the principal compenents, as well as a method to reduce the number of eigenvectors.

In [21]: pc = principal_components(img)
Covariance.....done

In [22]: w = view(pc.cov)

Image covariance matrix

In the covariance matrix display, whiter values indicate strong positive covariance, darker values indicate strong negative covariance, and grey values indicate covariance near zero.

To reduce dimensionality using principal components, we can sort the eigenvalues in descending order and then retain enough eigenvalues (an corresponding eigenvectors) to capture a desired fraction of the total image variance. We then reduce the dimensionality of the image pixels by projecting them onto the remaining eigenvectors. We will choose to retain a minimum of 99.9% of the total image variance.

In [23]: pc_0999 = pc.reduce(fraction=0.999)

In [24]: # How many eigenvalues are left?

In [25]: len(pc_0999.eigenvalues)
32

In [26]: img_pc = pc_0999.transform(img)

In [27]: w = view(img_pc[:,:,:3], stretch_all=True)

First 3 principal components of the image

Now we’ll use a Gaussian maximum likelihood classifier (GMLC) for the reduced principal components to train and classify against the training data.

In [28]: classes = create_training_classes(img_pc, gt)

In [29]: gmlc = GaussianClassifier(classes)
  Omitting class   7 : only 26 samples present
  Omitting class   9 : only 20 samples present
Covariance.....done
Covariance.....done
Covariance.....done
---// snip //---
Covariance.....done

In [30]: clmap = gmlc.classify_image(img_pc)
Classifying image...done

In [31]: clmap_training = clmap * (gt != 0)

In [32]: w = view_indexed(clmap_training)

GMLC results for 32 principal components

In [33]: training_errors = clmap_training * (clmap_training != gt)

In [34]: w = view_indexed(training_errors)

GMLC errors for 32 principal components

Fischer Linear Discriminant

The Fischer Linear Discriminant (a.k.a., canonical discriminant) attempts to find a set of transformed axes that maximize the ratio of the average distance between classes to the average distance between samples within each class. [Richards1999] This is written as

where is the covariance of the class centers and is the weighted average of the covariances of each class. If is invertible, this becomes

This eigenvalue problem is solved by the linear_discriminant function, yielding C-1 eigenvalues, where C is the number of classes.

In [35]: classes = create_training_classes(img, gt)

In [36]: fld = linear_discriminant(classes)
Covariance.....done
Covariance.....done
---// snip //---
Covariance.....done

In [37]: len(fld.eigenvectors)
15

Let’s view the image projected onto the top 3 components of the transform:

In [38]: img_fld = fld.transform(img)

In [39]: w = view_indexed(img_fld[:, :, :3])

Image data projected onto the first 3 components of the Fischer Linear Discriminant

Now we’ll classify the data using this discriminant.

In [40]: classes.transform(fld.transform)

In [41]: gmlc = GaussianClassifier(classes)

In [42]: clmap = gmlc.classify_image(img_fld)
Classifying image...done

In [43]: clmap_training = clmap * (gt != 0)

In [44]: w = view_indexed(clmap_training)

In [45]: fld_errors = clmap_training * (clmap_training != gt)

In [46]: w = view_indexed(fld_errors)

GMLC results for Fischer Linear Discriminant

See also

orthogonalize:

Gram-Schmidt Orthogonalization

Miscellaneous Functions

Band Resampling

Comparing spectra measured with a particular sensor to spectra collected by a different sensor often requires resampling spectra to a common band discretization. Spectral bands of a single sensor may drift enough over time such that spectra collected by the same sensor at different dates requires resampling.

For resampling purposes, SPy treats a sensor as having Gaussian spectral response functions for each of its spectral bands. A source sensor band will contribute to any destination band where there is overlap between the FWHM of the response functions of the two bands. If there is an overlap, an integral is performed over the region of overlap assuming the source band data value is constant over its FWHM (since we do not know the true spectral load over the source band) and the destination band has a Gaussian response function. Any target bands that do not have an overlapping source band will produce NaN as the resampled band value. If FWHM information is not available for a sensor’s bands, each band’s FWHM is assumed to reach half the distance the its adjacent bands. Resampling results are better when source bands are at a higher spectral resolution than the destination bands.

To create BandResampler object, we can either pass it a BandResampler object for each sensor or a list of band centers and, optionally, a list of FWHM values. Once the BandResampler is created, we can call it with a source sensor spectrum and it will return the resampled spectrum.

In [47]: import spectral.io.aviris as aviris

In [48]: img1 = aviris.open('f970619t01p02_r02_sc04.a.rfl', 'f970619t01p02_r02.a.spc')

In [49]: bands2 = aviris.read_aviris_bands('92AV3C.spc')

In [50]: resample = BandResampler(img1.bands, bands2)

In [51]: x1 = img1[96, 304]

In [52]: x2 = resample(x1)

NDVI

The Normalized Difference Vegetation Index (NDVI) is an indicator of the presence of vegetation. The index is commonly defined as

where NIR is the reflectance in the near infrared (NIR) part of the spectrum and RED is the reflectance of the red band. For our sample image, if we take band 21 (607.0 nm) for our red band and band 43 (802.5 nm) for near infrared, we get the following NDVI image.

In [53]: vi = ndvi(img, 21, 43)

In [54]: w = view(vi)

Normalized Difference Vegetation Index

ndvi is a simple convenience function. You could just as easily calculate the vegetation index yourself like this:

In [55]: red = img.read_band(21)

In [56]: nir = img.read_band(43)

In [57]: vi = (nir - red) / (nir + red)

Spectral Angles

A spectral angle refers to the angle between to spectra in N-space. In the absence of covariance data, spectral angles can be used for classifying data against a set of reference spectra by selecting the reference spectrum with which the unknown spectrum has the smallest angle.

To classify our sample image to the ground truth classes using spectral angles, we must compute the spectral angles for each pixel with each training class mean. This is done by the spectral_angles function. Before calling the function, we must first create a CxB array of training class mean spectra, where C is the number of training classes and B is the number of spectral bands.

In [58]: import numpy as np

In [59]: classes = create_training_classes(img, gt, True)
Covariance.....done
Covariance.....done
---// snip //---
Covariance.....done

In [60]: means = np.zeros((len(classes), img.shape[2]), float)

In [61]: for (i, c) in enumerate(classes):
   ....:      means[i] = c.stats.mean
   ....: 

In the code above, the True parameter to create_training_classes forces calculation of statistics for each class. Next, we call spectral_angles, which returns an MxNxC array, where M and N are the number of rows and columns in the image and there are C spectral angle for each pixel. To select the class with the smallest angle, we call the numpy argmin function to select the index for the smallest angle corresponding to each pixel. The clmap + 1 is used in the display command because our class IDs start at 1 (not 0).

In [62]: angles = spectral_angles(img, means)

In [63]: clmap = np.argmin(angles, 2)

In [64]: view_indexed((clmap + 1) * (gt != 0))

Spectral angle classification for training set

References

[Richards1999]

Richards, J.A. & Jia, X. Remote Sensing Digital Image Analysis: An Introduction. (Springer: Berlin, 1999).