Source code for spectral_libraries.core.music

# -*- coding: utf-8 -*-
"""
| ----------------------------------------------------------------------------------------------------------------------
| Date                : June 2020
| Copyright           : © 2020 by Arthur Maenhout (Locus) and Ann Crabbé (KU Leuven)
| Email               : acrabbe.foss@gmail.com
| Acknowledgements    : Jeroen Degerickx (KU Leuven). Documentation: https://doi.org/10.3390/rs9060565.
|
| This file is part of the Spectral Libraries QGIS plugin and python package.
|
| This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public
| License as published by the Free Software Foundation, either version 3 of the License, or any later version.
|
| This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied
| warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
|
| You should have received a copy of the GNU General Public License (COPYING.txt). If not see www.gnu.org/licenses.
| ----------------------------------------------------------------------------------------------------------------------
"""
import numpy as np


[docs]class Music: """ MUSIC (multiple signal classification) is a library pruning technique that operates on a full image. This technique reduces the dimensionality of an image to kf eigenvectors, and prunes a library according to the smallest distances to this subspace. Citation: Iordache, M.D., Bioucas-Dias, J.M., Plaza, A., Somers, B., 2014, MUSIC-CSR: Hyperspectral unmixing via multiple signal classification and collaborative sparse regression, IEEE Transactions on Geoscience and Remote Sensing, 52, p. 4364–4382. """ def __init__(self): self.log = None
[docs] @staticmethod def flatten_image(image: np.array) -> np.array: """ Store the image as a line (2D array). :param image: 3D image (m bands x n rows x p columns) :return: 2D image (m bands x n.p spectra """ n_bands = image.shape[0] n_pixels = image.shape[1] * image.shape[2] return np.reshape(image, [n_bands, n_pixels], order='F')
[docs] @staticmethod def normalize_brightness(spectra: np.array) -> np.array: """ Brightness normalization of spectral library or 2D-image :param spectra: array with m bands x n spectra :return: array with normalized spectral """ spectra_mean = np.mean(spectra, axis=0) normalized_spectra = np.divide(spectra, spectra_mean) return normalized_spectra
[docs] def estimate_noise(self, image): """ Hyperspectral noise estimation, by assuming that the reflectance at a given band is well modelled by a linear regression on the remaining bands. :param image: 2D hyperspectral data set (m bands x n pixels) :return: the noise estimates for every pixel (m x n) """ if image.shape[0] < 2: self.log('MUSIC: Too few bands to estimate the noise') return -1 small = 1e-6 n_bands, n_pixels = image.shape noise = np.zeros((n_bands, n_pixels)) rr = np.dot(image, np.transpose(image)) rri = np.linalg.inv(rr+small*np.eye(n_bands)) self.log('MUSIC: Computing noise for each band') for i in range(0, n_bands): part1 = rri - np.outer(rri[:, i], rri[i, :]) / rri[i, i] part2 = rr[:, i] part2[i] = 0 beta = np.matmul(part1, part2) beta[i] = 0 noise[i, :] = np.subtract(image[i, :], np.matmul(np.transpose(beta), image)) # noise_correlation = np.diag(np.diag(np.divide(np.matmul(noise, np.transpose(noise)), n_pixels))) return noise
[docs] def hysime2(self, image, min_eig): """ Hysime2 algorithm :param image: 2D hyperspectral data set (m bands x n pixels) :param min_eig: minimum number of eigenvectors to be retained from the image to calculate MUSIC distances :return: signal subspace dimension (tuple), eigenvectors that span the signal subspace (in the columns) """ noise = self.estimate_noise(image) # mxn n_bands, n_pixels = image.shape x = image - noise eigenvalue, eigenvector = np.linalg.eig(np.dot(x, x.T) / n_pixels) # the sign of eigenvectors changes from MATLAB to python, # since eigenvectors do not have a sign, this should not be relevant self.log('MUSIC: Estimating the number of endmembers') # calculate % variance explained by eigenvectors eigenvalue = eigenvalue / np.sum(eigenvalue) # [::-1] is added for descend sorting eigenvalue_sort = np.sort(eigenvalue)[::-1] eigenvalue_index = np.argsort(eigenvalue)[::-1] # sort the eigenvectors based on the sort of the eigenvalues eigenvector = eigenvector[:, eigenvalue_index] v_threshold = 0.999 v_sum = 0 kf = 0 while v_sum < v_threshold: v_sum = v_sum + eigenvalue_sort[kf] kf += 1 self.log('MUSIC: The first {} elements explain 99.9 percent of the variation'.format(kf)) if kf < min_eig: kf = min_eig self.log("MUSIC: Selecting the first {} eigenvectors".format(kf)) return eigenvector[:, 0:kf]
[docs] def execute(self, library: np.array, image: np.array, library_size: int, min_eig: int = 15, log: callable = print) -> dict: """ This technique reduces the dimensionality of an image to kf eigenvectors, and select library spectra according to the smallest distances to this subspace. :param library: the original spectral library to be pruned (bands x spectra), reflectance, no bad bands :param image: image used for pruning (bands x rows x columns) reflectance, no bad bands :param library_size: pruned library size :param min_eig: minimum number of eigenvectors to be retained from the image to calculate MUSIC distances :param log: communicate messages (point to the print_log tab in the GUI; otherwise print to the console) :return: dictionary with each output metric as a numpy array in a key-value pair - music_indices: indices of the pruned library - eigenvectors: first kf eigenvectors of image subspace (needed for sparse unmixing) kf = estimated number of eigenvectors in the image estimated by Hysime algorithm - music_distances: MUSIC distances of the pruned library spectra to the image, sorted from low to high """ self.log = log if log else print # reshape image to 2D image = self.flatten_image(image) # brightness normalization of library and image library = self.normalize_brightness(library) image = self.normalize_brightness(image) # delete all zero-spectra in the image image = np.delete(image, np.where(np.sum(image, axis=0) == 0)[0], axis=1) # hysime2 + final size of pruned library eigenvectors = self.hysime2(image, min_eig) # calculate distances from library spectra to image subspace n_bands = library.shape[0] # Euclidean distances between spectra and the hyperplane formed by the retained eigenvectors p = np.identity(n_bands) - np.matmul(eigenvectors, eigenvectors.T) error = np.sqrt((np.matmul(p, library) ** 2).sum(axis=0)) / np.sqrt((library ** 2).sum(axis=0)) # sort these distances error_sorted = np.sort(error) error_sorted_indices = np.argsort(error) library_indices = error_sorted_indices[0:library_size] self.log('MUSIC: {} spectra selected'.format(len(library_indices))) return {'music_indices': library_indices, 'music_eigenvectors': eigenvectors, 'music_distances': error_sorted[0:library_size]}