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. Chaffee County, Colorado. (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 Imaging and End Member Matching

The human eye is sensitive to three primary colors, red, green, and blue, broadly centered around wavelengths 560, 530, and 420nm, with considerable overlap. In contrast, hyperspectral imaging cameras, such as the Jet Propulsion Laboratory’s airborne AVIRIS sensor, may record up to 224 separate wavelengths in narrow bands between 370 (near UV), and 2500nm (near-thermal infrared) with little overlap, thereby affording detailed spectroscopy of reflected light and identification of the chemical composition of the reflecting surface. Hyperspectral or Multispectral Imaging may also feature in some MRI, CT, and LIDAR applications, but airborne and space-based remote sensing are the best-known examples.

Figure (1) is a ground photograph of Coulter’s field targets. 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 in Chaffee County 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 areas, 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. These three spectrally pure endmembers were determined by mineralogical estimation of the selected surface rocks, their spectra measured in the lab 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 research 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 in all applications. Mathematically, given a hyperpectral image matrix $F_{\text{ LxN}}$ of $N$ pixels and $L$ spectral channels, we wish to solve

for a mixture matrix Mix whose column vectors each represent the spectrum of a linear combination of distinct endmembers, and where $P\le L$ is the number of spectrally distinct endmembers in the scene. All coefficients of Mix must be non-negative, and the norm of each column vector is usually required to be 1. If each column of R is also so normalized, then likewise will be the Abundance matrix

The parameter $V_d$is termed the virtual dimensionality. Determining an appropriate value for a particlar image is part of the problem. The channel length $L$ can be excessively large: up to 224 for AVIRIS airborne cameras and 2048 for modern MRI. Dimensionality reduction techniques include slective subsampling, Principle Component Analysis and, if additional noise properties are known, Minimum Noise Fraction. Eq. (1) is then reduced to smaller row dimension:

where $P=S_d$, the subspace dimension with $P\le L$ and usually $P\ll L$. $S_d$and $V_d$are often the same.

Earlier approaches include the N-Findr algorithm proposed by Winter in 1999[4]. 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.

Further noise reduction is afforded by the nature of the SimFit algorithm itself: the simplex vertices are defined by the orientation and location of its faces, which are determined from a weighted average of pixels nearby to each face. The contributions from each member in the set of adjacent pixels are weighted by distance and summed, which tends to smooth out individual pixel noise.

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 solution 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 Initial Synthetic Image Tests

Brief descriptions of Icarus Image automated tests are listed in Appendix A. The most recent (as of 27 June 2024) tests the final SimFit Mix matrix against the original Synthetic Image’s Mix matrix 6 in the presence of $\pm 50\%$ added gaussian noise, or 100% positive gaussian noise, with variance $\sigma =10$. Nearly all components of final SimFit Mix matrix agree with those of the initial clean Synthetic Mix matrix – the clean mix matrix before noise was added – within 5% of maximum value.

This is a satisfying result. These initial synthetic images are completely random: they have no spatial structure and no spatial noise filtering was done.

The final SimFit Abundance matrix contains negative components corresponding to noisy pixels that lie outside the final fitted simplex. The 12% of Abundance columns that contain such negative components are set to zero. Although there are extremes near unity, nearly all the remaining non-negative SimFitAbundance coefficients lie within 1% of the original (noisy) Synthetic Abundance, and are representable by positive combinations of SimFit Mix matrix columns. It remains to decompose the final SimFit mixture matrix columns into their basis set components.

To this eventual end, as an internal consistency check this testcase also shows how SimFitMix = SynMix*X may be solved for decomposition matrix X via SVD. The decomposition is diagonal-dominant and accurate to within a few tenths of a percent in this test case where the original Synthetic Mix matrix and its dimensions are known. Beyond self-consistency and known original Mix Matrix, future tests will replace SVD with a constrained least-squares solution to satisfy real-world requirements for non-negative mixing coefficients. Analysis of the resulting residual matrices should be interesting.

These encouraging results suggest our next test suite:

1. Replace SVD placeholder with Constrained Least Squares.
2. Allow SimFit virtual_dimensions larger or smaller than the number of Synthetic endmembers.
3. Relax current testcase requirement that SimFit subspace_dimension = virtual_dimension. (This is a testcase limitation; IcarusImage itself allows subspace_dimension $\ge$ virtual_dimension.)
4. In analysis of converged SimFit Mix Matrix, allow test endmember basis sets not necessarily the same as the original Synthetic:
   • Test basis set smaller than the Synthetic basis, cannot span the actual spectral space.
   • Test basis larger than the Synthetic basis, assumes endmembers not in evidence.
   • Test basis of right size, but with one or more endmembers different from those in the synthetic.

   Such Constrained Least Squares tests should inform future endmember unmixing using Cross-Correlation matching, Constrained Energy Minimization, or Spectral Angle Mapping.

5. When a test SimFit Mix matrix with endmembers not in the Synthetic Image is projected into a noisy Synthetic Image’s PCA space, what to it’s PCA component coefficients look like regarding the Synthetic Image’s singular vectors and values?
6. How do these all behave as Signal to Noise Ratio drops beneath 1?
7. Add simple spatial structure to Synthetic Images, such as concentric rings, so one might visualize noise rejection as simplex optimization progresses.

Synthetic and phantom images incorporating real-world chemical surface or LIDAR reflectance, NMR, or X-Ray absorbtion endmember spectra, and sensor-specific noise models, might also be of interest.

6 Next Steps

The Simplex Fitting and optimization discussed here was an open problem when we first initiated this study. The N-Findr and Shrink-Wrap algorthms were introduced in 1999, Vertex Component Analyis in 2005. Upon these SimFit is an improvement. But we are not terribly unique, and SimFit or similar algorithms may well have been implemented by others in the interum. A thorough literature search, and requests for comments, 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. Explore implementing compatibility with data generated by NASA’s EMIT and PACE space-based sensors.
2. If there is interest, explore implementing compatibility with ENVI or other commercial products.
3. Compare VCA and SimFit with N-Finder, Shrink-Wrap, and more recent solutions when found.
4. 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 possible down-sampling 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 how one might optimize their selection?
5. 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 will require some effort, and might be unnecessary depending upon potential collaboration.
6. Explore possible applications to MRI, multi-energy CT, and multispectral LIDAR.

Improved dimensionality and noise[3] estimations should be considered throughout.

A Synthetic Image Tests

Testcases automated via the Qt Test module:

1. Pure Legendre Endmembers Tests with virtual_dim = 8 pure endmember image pixels that the generated VCA and Simfit Mixture and Abundance matrices are identical with their respective Synthetic image originals, within roundoff. These tests are trivial, but must be checked.
   • virtual_dim: 8
   • subspace_dim: 10
   • Nchannels: 256
   • NpurePixels: 8
   • PercentNoise: 0
   • Endmember_t: legendre
   • DimensionalityReduction_t: correlation

   This testcase is significant as it demonstrates expected Legendre polynomial orthogonality, that Vca algorithm behaves as expected with sufficient pure endmember image pixels, and that running SimFit following the Vca analysis does not change this result.

2. Impure Legendre Tests, with no pure endmember image pixels, that
   1. Original image R can be regenerated from Reduced Image X via Correlation Component back projection.
   2. Synthetic Mix matrix column inner products have expected Legendre polynomial orthogonality.
   3. SimFit Mix matrix can be regenerated from Legendre polynomials and the component matrix C of inner products of the original SimFit Mix matrix with Legendre polynomials, which is possible since the original Synthetic Mix matrix is composed of sums of two Legendre polynomials.
   4. The SimFit Mix component matrix C for the above test may be generated both directly from the inner products of SimFit Mix matrix columns with Legendre polynomials, AND by SVD solution of SimFitMix = PolyVal * C, where PolyVal is a 256x8 matrix for Legendre polynomials sampled on [-1, 1]

   These tests are relatively trivial, but must be checked.

   • virtual_dim: 8
   • subspace_dim: 10
   • Nchannels: 256
   • NpurePixels: 0
   • PercentNoise: 0
   • Endmember_t: legendre
   • DimensionalityReduction_t: correlation

   This testcase is significant as it demonstrates expected final SimFit Mixture matrix unmixing properties against endmember polynomials.

3. Impure Unitvec Tests final SimFit Mix and Abundance matrices against their Synthetic Image originals. Demonstrates the final SimFitMix is is the same as the original SynMix when the original Synthetic Image Mixture matrix is composed of orthogonal unit vectors. In this case each SynMix column $j=0,1,...7$ is unit vector of length 8 of values $\delta \left(j\right)$
   • virtual_dim: 8
   • subspace_dim: 8
   • Nchannels: 8
   • NpurePixels: 0
   • PercentNoise: 0
   • Endmember_t: unitvec
   • DimensionalityReduction_t: none

   Such unit vectors are orthogonal for any length, we choose Nchannels=8 for this test case.

   This test is significant as it demonstrates the final SimFit Mixture, Abundance, and Image matrices are identical within round-off to the Synthetic image originals, for noise-free synthetic images generated from unit vector endmembers, with no pure endmember image pixels. Here Vca alone fails as there are insufficient pure endmember image pixels.

4. Impure Synthetic Mix Tests final SimFit Mix and Abundance matrices against their Synthetic Image originals. Demonstrates the SimFitMix can be transformed back to the original SynMix when one knows the latter, which we do for synthetic tests. In this case each SynMix column $j=0,1,...7$ is a sum of two Legendre polynomials $L_0\left(x\right)+L_{j+1}\left(x\right)$ where $L_0\left(x\right)$ provides an offset to make all samples positive, and $x$ is sampled just 8 times on [-1,1]
   • virtual_dim: 8
   • subspace_dim: 8
   • Nchannels: 8
   • NpurePixels: 0
   • PercentNoise: 0
   • Endmember_t: legendre
   • DimensionalityReduction_t: none

   With so few channels our inner product integral returns true Legendre orthogonality only between the linear polynomials $L_0$ and $L_1$. But we can still use Eigen’s Svd.solve() to accurately compute SimFitMix unmixing against the original Synthetic image Mix matrix.

   This test is significant as it demonstrates the consistency of final SimFitMix unmixing for noise-free synthetic images with no pure endmember image pixels, a common case where Vca alone fails.

5. Impure Synthetic Mix with Noise Tests final SimFit Mix and Abundance matrices against their Synthetic Image originals that contain additive white gaussian noise. Demonstrates the SimFitMix can be transformed back to the original noisy synthetic SynMix when the problem is not underspecified, that is, when the number of channels is equal to the size of the endmember basis set. In this case each SynMix column $j=0,1,...7$ is a sum of Legendre polynomials $L_0\left(x\right)+L_{j+1}\left(x\right)$ where x is sampled just 8 times on [-1,1]
   • Nendmembers: 8
   • Nchannels: 8
   • virtual_dim: 8
   • subspace_dim: 8
   • NpurePixels: 0
   • PercentNoise: 5
   • Noise Sigma: 10
   • Endmember_t: legendre
   • DimensionalityReduction_t: none

   With so few channels our inner product integral returns true Legendre orthogonality only between. the linear polynomials $L_0$ and $L_1$. But with Nendmembers = Nchannels we can use Eigen’s Svd.solve() to accurately compute SimFitMix unmixing against the original Synthetic image Mix matrix. In this case the resulting SimFit Abundance matrix is the same as the original noisy Synthetic Abundance, and we can set to zero all Abundance columns which contain one or more negative coefficients as these represent noisy pixels that lie outside the fitted simplex.

   This test is significant as it demonstrates the consistency of this SimFitMix unmixing for noisy synthetic images with no pure endmember image pixels, a common case where Vca alone fails, and rejection of those noisy pixels that cannot be represented as linear combinations of the origian synthetic basis set. Many remaining pixels will still have noise not present in the original Synthetic Image Clean Abundance matrix, but these remaining pixels lie inside the fitted simplex and are representable by the original basis.

6. Impure Synthetic Mix with Noise and PCA/Correlation dimensionality reduction Tests final SimFit Mix and Abundance matrices against their Synthetic Image originals that contain additive white gaussian noise (awgn). Demonstrates the SimFitMix can be transformed back to the original noisy synthetic SynMix when the number of channels is large, 64 in this example, using pca via correlation to reduce the dimensionality to the known dimensionality of the synthetic image: 8 in this example. In this case each SynMix column $j=0,1,...7$ is a sum of Legendre polynomials $L_0\left(x\right)+L_{j+1}\left(x\right)$, where x is sampled 64 times on [-1,1]
   • Npixels: 2097152 (2048x1024)
   • Nendmembers: 8
   • Nchannels: 64
   • virtual_dim: 8
   • subspace_dim: 8
   • NpurePixels: 0
   • PercentNoise: 100
   • Noise Sigma: 10
   • Final alpha: 2.94e-3
   • Endmember_t: legendre
   • DimensionalityReduction_t: correlation

   With 64 channels our inner product integral returns expected Legendre orthogonality for all endmembers within 0.05%.

   In this testcase we tuned alpha_final to roughly zero the final SimFitMix.minCoeff(): min,max = (-6e-6, 2.000e+00). Larger alpha results in a larger more positive minimum, while smaller values yield a negative minimum. The root search to zero minCoeff() is applied algorithmically.

   In presence of 100% added gaussian noise nearly all components of the final SimFitMix agree with those of the initial clean Synthetic Mix matrix within 5% – the clean Synthetic mix *before* noise was added. Smaller alpha values include more pixels within the simplex, hence fewer with negative mixing coefficients.

   Although the extremes are near unity, nearly all non-negative SimFitAbundance coefficients lie within 1% of the original (noisy) Synthetic Abundance.

   As a consistency check, with pca and virtual_dim small, we use Eigen’s Svd.solve() to accurately decompose SimFitMix component matrix X with respect to the original Synthetic Mix in the reduced pca subspace. The component matrix X is then applied to the SimFit Abundance Matrix. In this case the resulting SimFit Abundance matrix product is the same as the original noisy Synthetic Abundance. We can again set to zero all Abundance columns which contain one or more negative coefficients, as these represent noisy pixels that lie outside the fitted simplex. In this test case 254643 negative coefficient noise pixels were removed from SimFitAbundance, 12% of the 2097152 pixels total. In practice coefficients within a small negative tolerance e.g. -1e-5 might be allowed.

   This test is significant as it demonstrates the internal consistency of this SimFitMix unmixing for noisy synthetic images with no pure endmember image pixels, a common case where Vca alone fails, and rejection of those noisy pixels that cannot be represented as linear combinations of the original synthetic basis set. Many remaining pixels will still have noise not present in the original Synthetic Image Clean Abundance matrix, but these remaining pixels lie inside the fitted simplex and are representable by the original basis.

   This test also demonstrates that for this sort of Synthetic, the number of channels in the Synthetic Image may greatly exceed VCA’s subspace_dimension, which nonetheless must still equal its virtual_dimension, for this test.

   In the special case where we wish to recover the original Image by back-projecting the PCA Reduced image, virtual_dim must equal Nchannels so that the full set of covariance eigenvectors is retained. However, such large virtual and subspace dimensions might introduce errors into the SimFit Mix matrix, and are unnecessary for the SimFit analysis. Here we used virtual_dim=subspace_dim=8, and Nchannels=64.

References