Figure 1: Ruby Mountain a.k.a. “West Red” as seen looking east from Enterprise Peak a.k.a. Larson’s Mountain. “East Red” (Mount Blaurock) in background, upper right. Garfield Peak is the craggy highpoint to the left of Ruby. (Photograph by author, August 2008)

Preface and Acknowledgement. Icarus Image and my own studies into hyperspectral image processing were inspired by David Coulter’s Remote Sensing Analysis of Alteration Mineralogy Associated with Natural Acid Drainage in the Grizzly Peak Caldera, Sawatch Range, Colorado, Colorado School of Mines doctoral dissertation, 2006. Apart from the obligatory Cuprite cameo (16), all aerial photography and imagery used herein were taken directly from Coulter’s dissertation, which represents a substantial effort on Coulter’s part organizing his acid-drainage study, conducting the ground mineralogical survey, and coordinating the remote-sensing overflights by the Jet Propulsion Laboratory’s AVIRIS aircraft and SpectIR Corporation’s HST instrument. My own desk-efforts in no way compare with David Coulter’s.

Contents

1 Introducing Hyperspectral End Member Matching

Figure (1) is a ground photograph of the target of this exercise. Ruby Mountain (aka Red Mountain or “West Red”) and its “East Red” companion lie in the southwest part of Colorado’s Grizzly Caldera, an extensive volcanic feature fifteen miles southeast of Aspen. Ruby Mountain and mine are readily found on Google Earth.¹ Photo was taken from Enterprise Peak (Larson’s Mountain), on the Sawatch Range separating Grizzly Caldera from Taylor Park. Foreground drainage is the southern end of Lincoln Gulch, which drains north (left) to the Roaring Fork on the western approach to Independence Pass. The distinctive orange color is characteristic of Hematite (a.k.a. “rust”), the primary natural oxidation state of iron. But there are contributions from Goethite and Jarosite as well, which form when iron weathers under more acidic conditions. Coulter sought to map, by remote sensing imagery, the geographic distribution of these three iron oxides in the target area, and thereby characterize naturally occuring acid seeps.

Figure (2) represents a simple example of supervised end-member classification. In this case the supervisor (David Coulter) painstakingly took “ground-truth” samples from select locations in the target area, and identified Hematite, Goethite, and Jarosite as the end-members of interest. Spectrally pure endmembers were determined by mineralogical estimation of the selected surface rocks, and their spectra measured and recorded for later matching against many (possibly millions) such pixels in the hyperspectral images.

The target pixel’s spectrum is shown in red, the spectrum resulting from the best fits of the three end members by least-squares fitting and by our implementation of the Optimized Cross-Correlation Matching (OCCM) algorithm used by Coulter are shown in blue and green respectively, with vertical offsets added for clarity. The lowest-most curve illustrates goodness-of-fit (Pearson’s $R=0.998$). Matching was done over two wavelength bands: 370-1100 nm, and 2000-2440 nm. The intervening 1100-2000 nm range is characterized by highly variable water absorbtion, and is frequently excluded from this sort of study. (Hence the figure’s flat plateau in this region.)

2 Icarus Image VCA: Vertex Component Analysis with Simplex Refinement

In contrast to OCCM, Icarus Image is part of an ongoing project to investigate algorithms for unsupervised hyperspectral endmember classification, that is, how to identify and map chemical endmember spectra in an image in absence of local ground-truth sampling. This problem is more difficult, for without ground truth a wider selection of endmember spectra must be applied to match the spectra of each pixel in the image. Restricting the endmember set to a the (relatively) few compounds appropriate to a particular image region becomes an interesting problem.

A common approach, which we follow here, is to first estimate the number M of spectrally distinct endmembers in the area of interest – itself a good problem – then attempt to find a minimal-volume M-dimensional simplex in spectral space that best contains all (or nearly all) the individual spectra from each image pixel. Under the linear-mixing assumption the spectra of each pixel contained within the simplex may then be expressed through a (relatively) low-dimension mixing matrix as a linear combination of the spectra represented by the vertices of the simplex. Linear mixing is generally valid for the remote sensing reflectance spectroscopy considered here, but might not be for all cases, for instance in functional MRI.

Earlier approaches include the N-Findr algorithm proposed by Winter in 1999[3]. N-Findr finds a maximum-volume simplex whose vertices lie within the set of image pixels. A related approach is Vertex Component Analyis (VCA)[2], which we incorporate here as an intermediate step.

Related to our final simplex optimization (SimFit) is the Shrink-Wrap approach discussed by Fuhrmann[1]. But although our goals are similar, the two are not the same. In particular, Shrink-Wrap reportedly can yield simplices that completely contain all the pixel spectra in the area of interest, which is not necessarily what one desires. First there is sensor noise. Even a few mis-read pixels may result in an unsuitably large simplex with vertices in spectral space lying far outside the chemically-realizable spectral set of interest. Second is the possibility of a relatively small number of pixels from endmember compounds that either are not of interest and could expand the simplex in directions not appropriate to the area at large, or which are of interest but whose small number warrant separate independent analysis and classification.

Icarus Image is a research project whose principal algorithmic goal is to find a minimal-volume spectral simplex that contains the spectra of almost all the pixels in the area of interest, and can identify those pixels that lie far outside. “Almost” and “far” are quantifiable parameters of the problem. As with nearly all M-dimensional optimizers, our minimal-volume solution simplices will not be unique, but for a given dataset each should yield nearly the same final chemical end-member mixing matrix.

To this end Icarus Image’s graphical interface serves as a flexible test platform, enabling rapid changes in program option editing and display. At this point (November 2023) the GUI is completely standalone, with ability to export image and spectral signature files to Ball Aerospace’ Opticks GIS.

Programmatically, our VCA (Vertex Component Analysis) and SimFit (simplex optimization and refinement) functionalities are implemented in a shared code library, for which there are separate GUI and command-line drivers. The latter allows rapid and scripted execution of simple test cases for development, and batch processing when desired. The Qt GUI introduces flexibility into the command line’s otherwise linear control flow, and enables convenient display of the input and output images. Briefly

1. An input hyperspectral image is selected. Format may be either ENVI Standard or GeoTiff.
2. Optionally, the coordinates of a desired sub-image of interest may be given either on the command line, in a configuration file, or set in the GUI’s Options Editor.
3. Optionally, a suitably georeferenced rgb or greyscale background aerial photo covering the same areas of interest may also be loaded. Rgb bands from the hyperspectral image or its derived PCA, VCA, or SimFit images may then be be overlayed and alpha-blended over the background image.
4. Dimensionality Reduction. We presently incorporate three traditional eigenthresholding based methods: Harsanyi-Farrand-Chang (HFC), noise-whitened HFC (NWHFC), and noise subspace projection (NSP) to determine the first M PCA bands most likely to represent the major hyperspectral endmembers. VCA-related Othogonal Subspace Projection (OSP), and Simplex Fitting-related condition-number based methods, are also under consideration. Dimensionalty-reduction possibilities are far from exhausted.

   
   In the case of lower-dimension multispectral images, such as from the 16-band World View 3 satellite, or by downsampling hyperspectral images to multispectral dimensions, first-step PCA analysis is still possible if not strictly necessary. (Such down-sampling capability is in progress).

5. Run Vertex Component Analyis (VCA) to get an initial spectral simplex with vertices associated with specific and identifiable hyperspectral image pixels. Classic VCA starts with an arbitrary direction in spectral space that is not orthogonal to any image pixel spectra, then locates the suitably-projected image pixel spectral with the greatest projection in this chosen spectral direction. This pixel is the first VCA component vertex.

   
   The second VCA vertex is detemined by choosing a second spectral dimension orthogonal to the spectral direction of the first pixel vertex, and finding the next projected pixel with greatest projection in this second spectral direction. A third spectral direction is choosen orthogonal to the spectral direction of the first two, and a third VCA vertex pixel found as that with greatest projected projection in this third direction.

   The remaining M-3 VCA vertices are similarly found by recursion.

   Vertex Component Analysis is rapid, and Really Neat – provided there neatly exist M image pixels that represent pure and distinct spectral (and chemical) endmembers. Here on planet earth such neatness is rarely obtained save in carefully constructed synthetic images. Most natural hyperspectral images will contain a large number of pixels whose spectral coordinates lie outside (usually far outside) the M-1 dimensional simplex defined by the M VCA coordinates.

   Our goal here is to find a Mixing Matrix that defines the spectrum of each image pixel in terms of the spectra of the simplex vertices. But this is only meaningful for image pixels whose spectra lie inside the spectral simplex whose vertices are defined by the M VCA vertices. So while VCA is intuitively appealing, it somewhat lacks in practical geophysical application. Therefore...

6. Run Simplex Fitting optimization and refinement, which expands the spectral coordinates of the VCA simplex vertices such that the resulting spectral simplex contains most of the image pixels, where “most” is characterized by an input tuning parameter $\alpha$. In this respect the SimFit algorithm used may be thought of as an improved Shrink-Wrap, skirting the shrink-wrap issues mentioned above.

   
   SimFit finds an optimal simplex with some pixels’ spectra lying on each side of each simplex face. By far most pixels will lie inside the final simplex when converged, and most of the (relatively few) outside pixels will not lie far outside, as illustrated in Figures (10) and (11). The remaining (presumably very few) far outliers, whether noise or otherwise, are another story. Their characterization is TBD.

Icarus Image screenshots are illustrated in Figures (3) through (9):

3 Prelude to Spectroscopy

2015 efforts included refining the refinement algorithm, noting some limitations, and contemplating future modifications to avoid them. Some cpu-intensive portions were moved to custom CUDA kernels. GPU speedup was good enough to suggest other approaches were well within grasp. Year 2022 saw substantial speedups to both CPU and GPU portions of the code, and incorporation of an independent gradient optimizer for confirmation and comparison.

Spectral matching has progressed with introduction of Ball Aerospace’s Opticks GIS. Experience with the earlier OCCM supervised matching algorithms afforded some insight into what an automated matching program might entail; I was relieved to find a robust open-source GIS that will largely fill our needs.

Opticks is written with specific focus on multi- and hyper- spectral image analysis, with underlying image data structure and file formats tailored to that end. Its Spectral Plugin and Library Builder will likely be adequate as-is, after convolution of JPL/USGS spectral libraries to the resolution and band shapes of AVIRIS and other sensors. Some early results:

4 Picture Show

We’ve done nearly all our development work referencing the small oxide-rich Ruby Mountain area with AVIRIS data illustrated above. This section includes a few other datasets from the Grizzly Caldera – including some from the SpectIR HST sensor – and the public AVIRIS image from the Cuprite district in Nevada. AVIRIS and HST images both have pixel size 2.5m x 2.5m. This is not analysis, which must await spectral matching capability. We merely illustrate SimFit robustness, and give some flavor for execution speed.

To this end all images in this section were processed using 16 PCA channels and 16 Virtual Dimensions, resulting in 16x16 Mixture Matrices. There were 122 AVIRIS bands in the wavelength intervals of interst. The workstation has a 16-core 3.4 GHz AMD-5950X processor, 128 GB ECC ram, and Nvidia Quadro RTX 4000 (Turing) 8GB graphics. We build both 32-bit and 64-bit floating point (fp32 and fp64) executables. Only the 32-bit version uses the GPU, affording approximately 4x speedup over fp64 for the larger images.

5 Next Steps

The Simplex Fitting and optimization discussed here was an open problem when we initiated this study fifteen years ago. The N-Findr and Shrink-Wrap algorthms were introduced in 1999, Vertex Component Analyis in 2005. Upon these SimFit is an improvement. But I am not terribly unique, and SimFit or similar algorithms may well have been implemented by others in the interum. A thorough literature search, and comment requests, are much in order before the Icarus Image project proceeds much further.

But in the event Icarus Image does offer worthwhile additions, the following might be pursued:

1. Add automated unit and synthetic image tests. Our Synthetic Image subproject showed positive results some years ago; it needs expansion and incorporation into an automated test suite.
2. Explore implementing compatibility with data generated by NASA’s EMIT and PACE space-based sensors.
3. If there is interest, explore implementing compatibility with ENVI or other commercial products.
4. Compare VCA and SimFit with N-Finder, Shrink-Wrap, and more recent solutions as they turn up.
5. Finish capability to down-sample hyperspectral (224 channel AVIRIS and HST) images to lower-dimension Multispectral (e.g. 16 channel World View 3, more for newer sensors) as alternative to PCA dimensionality reduction. Explore possibile downsampling to arbitrary Multispectral dimensions with band widths and spacing adjustable on a per-image basis. If such band adjustments can improve the resolution of particlar features, can one visualize a cost-function to optimize their selection?
6. Spectroscopy: write USGS / JPL spectral database interfaces to Opticks’ Spectral spectral-matching plugin and Library Builder. Government and commercial hyperspectral analysts already have production-quality spectral databases. Reimplementing our own might require some effort, but proof-of-concept for synthetic images mightn’t be unreasonable.
7. Explore possible applications to MRI.

Improved dimensionality estimation should be considered throughout.

References