# Source code for pysptools.eea.eea_int

#
#------------------------------------------------------------------------------
# Copyright (c) 2013-2014, Christian Therien
#
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#
# Unless required by applicable law or agreed to in writing, software
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
#------------------------------------------------------------------------------
#
# eea_int.py - This file is part of the PySptools package.
#

"""
PPI, NFINDR, ATGP, FIPPI classes
"""

import numpy as np
from . import eea
from .inval import *
from .docstring import *

def _normalize(M):
"""
Normalizes M to be in range [0, 1].

Parameters:
M: numpy array
1D, 2D or 3D data.

Returns: numpy array
Normalized data.
"""
minVal = np.min(M)
maxVal = np.max(M)

Mn = M - minVal;

if maxVal == minVal:
return np.zeros(M.shape);
else:
return Mn / (maxVal-minVal)

def _plot_end_members(path, E, utype, is_normalized, axes, suffix):
""" Plot a endmembers graph using matplotlib """
import os.path as osp
import matplotlib.pyplot as plt
if axes == None:
axes = {}
axes['wavelength'] = [x+1 for x in range(E.shape[1])]
axes['x'] = 'Wavelength'
axes['y'] = 'Brightness'
else:
if not('wavelength' in axes) or axes['wavelength'] == None:
axes['wavelength'] = [x+1 for x in range(E.shape[1])]
if not('x' in axes) or axes['x'] == None:
axes['x'] = 'Wavelength'
if not('y' in axes) or axes['y'] == None:
axes['y'] = 'Brightness'

plt.ioff()
plt.xlabel(axes['x'])
if is_normalized == True:
plt.ylabel(axes['y']+' - normalized')
else:
plt.ylabel(axes['y'])
plt.title('Spectral Profile')
plt.grid(True)
n_graph = 1
legend = []
for i in range(E.shape[0]):
plt.plot(axes['wavelength'], E[i])
legend.append('EM{0}'.format(str(i+1)))
if (i+1) % 5 == 0 :
plt.legend(legend, loc='upper left', framealpha=0.5)
legend = []
if suffix == None:
fout = osp.join(path, 'emembers_{0}__{1}.png'.format(utype, n_graph))
else:
fout = osp.join(path, 'emembers_{0}__{1}_{2}.png'.format(utype, n_graph, suffix))
try:
plt.savefig(fout)
except IOError:
raise IOError('in _plot_end_members, no such file or directory: {0}'.format(path))
n_graph += 1
plt.clf()
plt.xlabel(axes['x'])
if is_normalized == True:
plt.ylabel(axes['y']+' - normalized')
else:
plt.ylabel(axes['y'])
plt.title('Spectral Profile')
plt.grid(True)
if E.shape[0] % 5 != 0:
plt.legend(legend, loc='upper left', framealpha=0.5)
if suffix == None:
fout = osp.join(path, 'emembers_{0}__{1}.png'.format(utype, n_graph))
else:
fout = osp.join(path, 'emembers_{0}__{1}_{2}.png'.format(utype, n_graph, suffix))
try:
plt.savefig(fout)
except IOError:
raise IOError('in _plot_end_members, no such file or directory: {0}'.format(path))
plt.close()

def _display_end_members(U, utype, is_normalized, axes, suffix):
""" Display endmembers using matplotlib to the IPython Notebook. """
import matplotlib.pyplot as plt
if axes == None:
axes = {}
axes['wavelength'] = [x+1 for x in range(U.shape[1])]
axes['x'] = 'Wavelength'
axes['y'] = 'Brightness'
else:
if not('wavelength' in axes) or axes['wavelength'] == None:
axes['wavelength'] = [x+1 for x in range(U.shape[1])]
if not('x' in axes) or axes['x'] == None:
axes['x'] = 'Wavelength'
if not('y' in axes) or axes['y'] == None:
axes['y'] = 'Brightness'

plt.xlabel(axes['x'])
if is_normalized == True:
plt.ylabel(axes['y']+' - normalized')
else:
plt.ylabel(axes['y'])
n_graph = 1
legend = []
for i in range(U.shape[0]):
plt.plot(axes['wavelength'], U[i])
legend.append('EM{0}'.format(str(i+1)))
if (i+1) % 5 == 0 :
plt.legend(legend, loc='upper left', framealpha=0.5)
legend = []
plt.xlabel(axes['x'])
if is_normalized == True:
plt.ylabel(axes['y']+' - normalized')
else:
plt.ylabel(axes['y'])
if suffix == None:
plt.title('Spectral Profile {0} - {1}'.format(n_graph, utype))
else:
plt.title('Spectral Profile {0} - {1} - {2}'.format(n_graph, utype, suffix))
plt.grid(True)
plt.show()
plt.clf()
n_graph += 1
if U.shape[0] % 5 != 0:
plt.legend(legend, loc='upper left', framealpha=0.5)
plt.xlabel(axes['x'])
if is_normalized == True:
plt.ylabel(axes['y']+' - normalized')
else:
plt.ylabel(axes['y'])
if suffix == None:
plt.title('Spectral Profile {0} - {1}'.format(n_graph, utype))
else:
plt.title('Spectral Profile {0} - {1} - {2}'.format(n_graph, utype, suffix))
plt.grid(True)
plt.show()
plt.close()

def _document(cls):
import sys
if sys.version_info[0] == 2:
cls.plot.__func__.__doc__ = plot_docstring
cls.display.__func__.__doc__ = display_docstring
if sys.version_info[0] == 3:
cls.plot.__doc__ = plot_docstring
cls.display.__doc__ = display_docstring

cmp = np.ndarray((n, vec.shape[1]), dtype=np.float)
i = 0
cmp[i] = vec[j]
i += 1
return cmp

[docs]class PPI(object): """ Performs the pixel purity index algorithm for endmembers finding. """ def __init__(self): self.E = None self.w = None self.idx = None self.idx3D = None self.is_normalized = False @ExtractInputValidation1('PPI') def extract(self, M, q, numSkewers=10000, normalize=False, mask=None): """ Extract the endmembers. Parameters: M: numpy array A HSI cube (m x n x p). q: int Number of endmembers to find. numSkewers: int [default 10000] Number of "skewer" vectors to project data onto. In general, recommendation from the literature is 10000 skewers. mask: numpy array [default None] A binary mask, when *True* the corresponding signal is part of the endmembers search. Returns: numpy array Recovered endmembers (N x p). """ if normalize == True: M = _normalize(M) self.is_normalized = True h, w, numBands = M.shape self.h, self.w, self.numBands = M.shape self.q = q M = np.reshape(M, (self.w*h, M.shape[2])) if isinstance(mask, np.ndarray): m = np.reshape(mask, (self.w*h)) cM = _compress(M, m) else: cM = M self.E, self.idx = eea.PPI(cM, q, numSkewers) self.idx3D = [(i % self.w, i // self.w) for i in self.idx] return self.E def __str__(self): return 'pysptools.eea.eea_int.PPI object, hcube: {0}x{1}x{2}, n endmembers: {3}'.format(self.h, self.w, self.numBands, self.q)
[docs] def get_idx(self): """ Returns: numpy array Array of indices into the HSI cube corresponding to the induced endmembers """ return self.idx3D
@PlotInputValidation('PPI') def plot(self, path, axes=None, suffix=None): _plot_end_members(path, self.E, 'PPI', self.is_normalized, axes=axes, suffix=suffix) @DisplayInputValidation('PPI') def display(self, axes=None, suffix=None): _display_end_members(self.E, 'PPI', self.is_normalized, axes=axes, suffix=suffix)
_document(PPI)
[docs]class NFINDR(object): """ N-FINDR endmembers induction algorithm. """ def __init__(self): self.E = None self.Et = None self.w = None self.idx = None self.it = None self.idx3D = None self.is_normalized = False @ExtractInputValidation2('NFINDR') def extract(self, M, q, transform=None, maxit=None, normalize=False, ATGP_init=False, mask=None): """ Extract the endmembers. Parameters: M: numpy array A HSI cube (m x n x p). q: int The number of endmembers to be induced. transform: numpy array [default None] The transformed 'M' cube by MNF (m x n x components). In this case the number of components must == q-1. If None, the built-in call to PCA is used to transform M in q-1 components. maxit: int [default None] The maximum number of iterations. Default is 3*p. normalize: boolean [default False] If True, M is normalized before doing the endmembers induction. ATGP_init: boolean [default False] Use ATGP to generate the first endmembers set instead of a random selection. mask: numpy array [default None] A binary mask, when *True* the corresponding signal is part of the endmembers search. Returns: numpy array Set of induced endmembers (N x p). References: Winter, M. E., "N-FINDR: an algorithm for fast autonomous spectral end-member determination in hyperspectral data", presented at the Imaging Spectrometry V, Denver, CO, USA, 1999, vol. 3753, pgs. 266-275. Note: The division by (factorial(p-1)) is an invariant for this algorithm, for this reason it is skipped. """ from . import nfindr if normalize == True: M = _normalize(M) self.is_normalized = True h, w, numBands = M.shape self.h, self.w, self.numBands = M.shape self.q = q M = np.reshape(M, (self.w*h, M.shape[2])) if transform != None: transform = np.reshape(transform, (self.w*h, transform.shape[2])) if isinstance(mask, np.ndarray): m = np.reshape(mask, (self.w*h)) cM = _compress(M, m) else: cM = M self.E, self.Et, self.idx, self.it = nfindr.NFINDR(cM, q, transform, maxit, ATGP_init) self.idx3D = [(i % self.w, i // self.w) for i in self.idx] return self.E def __str__(self): return 'pysptools.eea.eea_int.NFINDR object, hcube: {0}x{1}x{2}, n endmembers: {3}'.format(self.h, self.w, self.numBands, self.q)
[docs] def get_idx(self): """ Returns : numpy array Array of indices into the HSI cube corresponding to the induced endmembers """ return self.idx3D
[docs] def get_iterations(self): """ Returns : int The number of iterations. """ return self.it
def get_endmembers_transform(self): return self.Et @PlotInputValidation('NFINDR') def plot(self, path, axes=None, suffix=None): _plot_end_members(path, self.E, 'NFINDR', self.is_normalized, axes=axes, suffix=suffix) @DisplayInputValidation('NFINDR') def display(self, axes=None, suffix=None): _display_end_members(self.E, 'NFINDR', self.is_normalized, axes=axes, suffix=suffix)
_document(NFINDR)
[docs]class ATGP(object): """ Automatic target generation process endmembers induction algorithm. """ def __init__(self): self.E = None self.w = None self.idx = None self.idx3D = None self.is_normalized = False @ExtractInputValidation3('ATGP') def extract(self, M, q, normalize=False, mask=None): """ Extract the endmembers. Parameters: M: numpy array A HSI cube (m x n x p). q: int Number of endmembers to be induced (positive integer > 0). normalize: boolean [default False] Normalize M before unmixing. mask: numpy array [default None] A binary mask, if True the corresponding signal is part of the endmembers search. Returns: numpy array Set of induced endmembers (N x p). References: A. Plaza y C.-I. Chang, "Impact of Initialization on Design of Endmember Extraction Algorithms", Geoscience and Remote Sensing, IEEE Transactions on, vol. 44, no. 11, pgs. 3397-3407, 2006. """ if normalize == True: M = _normalize(M) self.is_normalized = True h, w, numBands = M.shape self.h, self.w, self.numBands = M.shape self.q = q M = np.reshape(M, (self.w*h, M.shape[2])) if isinstance(mask, np.ndarray): m = np.reshape(mask, (self.w*h)) cM = _compress(M, m) else: cM = M self.E, self.idx = eea.ATGP(cM, q) self.idx3D = [(i % self.w, i // self.w) for i in self.idx] return self.E def __str__(self): return 'pysptools.eea.eea_int.ATGP object, hcube: {0}x{1}x{2}, n endmembers: {3}'.format(self.h, self.w, self.numBands, self.q)
[docs] def get_idx(self): """ Returns: numpy array Array of indices into the HSI cube corresponding to the induced endmembers """ return self.idx3D
@PlotInputValidation('ATGP') def plot(self, path, axes=None, suffix=None): _plot_end_members(path, self.E, 'ATGP', self.is_normalized, axes=axes, suffix=suffix) @DisplayInputValidation('ATGP') def display(self, axes=None, suffix=None): _display_end_members(self.E, 'ATGP', self.is_normalized, axes=axes, suffix=suffix)
_document(ATGP)
[docs]class FIPPI(object): """ Fast Iterative Pixel Purity Index (FIPPI) endmembers induction algorithm. """ def __init__(self): self.E = None self.w = None self.idx = None self.idx3D = None self.is_normalized = False @ExtractInputValidation4('FIPPI') def extract(self, M, q=None, maxit=None, normalize=False, mask=None): """ Extract the endmembers. Parameters: M: numpy array A HSI cube (m x n x p). q: int [default None] Number of endmembers to be induced, if None use HfcVd to determine the number of endmembers to induce. maxit: int [default None] Maximum number of iterations. Default = 3*q. normalize: boolean [default False] Normalize M before unmixing. mask: numpy array [default None] A binary mask, when *True* the corresponding signal is part of the endmembers search. Returns: numpy array Set of induced endmembers (N x p). References: Chang, C.-I., "A fast iterative algorithm for implementation of pixel purity index", Geoscience and Remote Sensing Letters, IEEE, vol. 3, no. 1, pags. 63-67, 2006. """ if normalize == True: M = _normalize(M) self.is_normalized = True h, w, numBands = M.shape self.h, self.w, self.numBands = M.shape self.q = q M = np.reshape(M, (self.w*h, M.shape[2])) if isinstance(mask, np.ndarray): m = np.reshape(mask, (self.w*h)) cM = _compress(M, m) else: cM = M self.E, self.idx = eea.FIPPI(cM, q=q, maxit=maxit) self.idx3D = [(i % self.w, i // self.w) for i in self.idx] return self.E def __str__(self): return 'pysptools.eea.eea_int.FIPPI object, hcube: {0}x{1}x{2}, n endmembers: {3}'.format(self.h, self.w, self.numBands, self.q)
[docs] def get_idx(self): """ Returns: numpy array Array of indices into the HSI cube corresponding to the induced endmembers. """ return self.idx3D
@PlotInputValidation('FIPPI') def plot(self, path, axes=None, suffix=None): _plot_end_members(path, self.E, 'FIPPI', self.is_normalized, axes=axes, suffix=suffix) @DisplayInputValidation('FIPPI') def display(self, axes=None, suffix=None): _display_end_members(self.E, 'FIPPI', self.is_normalized, axes=axes, suffix=suffix)
_document(FIPPI)