Source code for pysptools.abundance_maps.amaps_int

#
#------------------------------------------------------------------------------
# Copyright (c) 2013-2014, Christian Therien
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#    http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
#------------------------------------------------------------------------------
#
# amaps_int.py - This file is part of the PySptools package.
#

"""
UCLS, NNLS, FCLS classes
"""

import os.path as osp
import numpy as np
from . import amaps
from .inval 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_abundance_map(amap, path, map_type, mask, interpolation, colorMap, columns, suffix):
    """ Plot an abundance map using matplotlib """

    def one_figure(amap, path, map_type, mask, interpolation, colorMap, columns, suffix):
        rows, remainder = divmod(amap.shape[2], columns)
        if remainder != 0: rows += 1
        n_amaps = amap.shape[2]
        f, axes = plt.subplots(rows, columns, figsize=(4*columns, 4*rows))
        f.set_dpi(128)
        #f.set_dpi(64)
        if suffix == None:
            f.suptitle('{0} Inversion'.format(map_type), fontsize=12)
        else:
            f.suptitle('{0} Inversion {1}'.format(map_type, suffix), fontsize=12)
        i = 1
        for r in range(rows):
            for c in range(columns):
                axes[r,c].axis('off')
        for r in range(rows):
            for c in range(columns):
                if i > n_amaps: continue
                m = amap[:,:,i-1]
                if isinstance(mask, np.ndarray):
                    m = m * mask
                ax = axes[r,c].imshow(m, interpolation=interpolation)
                ax.set_cmap(colorMap)
                axes[r,c].set_title('EM{}'.format(i), fontsize=10)
                i += 1
        if path != None:
            if suffix == None:
                fout = osp.join(path, '{0}.png'.format(map_type))
            else:
                fout = osp.join(path, '{0}_{1}.png'.format(map_type, suffix))
            try:
#                f.savefig(fout, dpi='figure')
                f.savefig(fout)
            except IOError:
                raise IOError('in abundance_map._plot_abundance_map, no such file or directory: {0}'.format(path))
        else:
            plt.show()
        plt.clf()

    def multi_figures(amap, path, map_type, mask, interpolation, colorMap, suffix):
        for i in range(amap.shape[2]):
            m = amap[:,:,i]

            if isinstance(mask, np.ndarray):
                m = m * mask

            img = plt.imshow(m, interpolation=interpolation)
            img.set_cmap(colorMap)
            plt.colorbar()
            if path != None:
                if suffix == None:
                    fout = osp.join(path, '{0}_{1}.png'.format(map_type, i+1))
                else:
                    fout = osp.join(path, '{0}_{1}_{2}.png'.format(map_type, i+1, suffix))
                try:
                    plt.savefig(fout)
                except IOError:
                    raise IOError('in abundance_map._plot_abundance_map, no such file or directory: {0}'.format(path))
            else:
                if suffix == None:
                    plt.title('{0} Inversion - EM{1}'.format(map_type, i+1))
                else:
                    plt.title('{0} Inversion - EM{1} - {2}'.format(map_type, i+1, suffix))
                plt.show()
            plt.clf()

    import matplotlib.pyplot as plt
    import warnings
    if path != None:
        plt.ioff()
    if columns != None and columns >= amap.shape[2]:
         warnings.warn('In abundance_map._plot_abundance_map, the number of abundances map to display is less or equal the number of columns')
    if columns != None and columns < amap.shape[2]:
        one_figure(amap, path, map_type, mask, interpolation, colorMap, columns, suffix)
    else:
        multi_figures(amap, path, map_type, mask, interpolation, colorMap, suffix)


_plot_docstring = """
        Plot the abundance maps.

        Parameters:
            path: `string`
              The path where to put the plot.

            mask: `numpy array [default None]`
              A binary mask, when *True* the selected pixel is displayed.

            interpolation: `string [default none]`
              A matplotlib interpolation method.

            colorMap: `string [default jet]`
              A matplotlib color map.

            columns: `int [default None]`
              Display all the images in one figure organized by
              columns.

            suffix: `string [default None]`
              Suffix to add to the file name.
        """


_display_docstring = """
        Display the abundance maps to a IPython Notebook.

        Parameters:
            mask: `numpy array [default None]`
                A binary mask, when *True* the selected pixel is displayed.

            interpolation: `string [default none]`
              A matplotlib interpolation method.

            colorMap: `string [default jet]`
              A matplotlib color map.

            columns: `int [default None]`
              Display all the images in one figure organized by
              columns.

            suffix: `string [default None]`
              Suffix to add to the title.
        """


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


def _compress(vec, mask):
    n = np.sum(mask)
    cmp = np.ndarray((n, vec.shape[1]), dtype=np.float)
    i = 0
    for j in range(mask.shape[0]):
        if mask[j] == 1:
            cmp[i] = vec[j]
            i += 1
    return cmp


def _expand(amap, mask, l, q):
    exp = np.zeros((l,q), dtype=np.float)
    i = 0
    for j in range(mask.shape[0]):
        if mask[j] == 1:
            exp[j] = amap[i]
            i += 1
    return exp


[docs]class UCLS(object): """ Performs unconstrained least squares abundance estimation. """ def __init__(self): self.amap = None self.m = None self.n = None self.q = None @MapInputValidation('UCLS') def map(self, M, U, normalize=False, mask=None): """ Performs unconstrained least squares abundance estimation on the HSI cube M using the spectral library U. Parameters: M: `numpy array` A HSI cube (m x n x p). U: `numpy array` A spectral library of endmembers (q x p). normalize: `boolean [default False]` If True, M and U are normalized before doing the spectra mapping. mask: `numpy array [default None]` A binary mask, when *True* the selected pixel is unmixed. Returns: `numpy array` An abundance maps (m x n x q). """ h,w,numBands = M.shape if normalize == True: M = _normalize(M) U = _normalize(U) Mr = np.reshape(M, (w*h, numBands)) if isinstance(mask, np.ndarray): m = np.reshape(mask, (w*h)) cMr = _compress(Mr, m) c_amap2D = amaps.UCLS(cMr, U) amap2D = _expand(c_amap2D, m, w*h, U.shape[0]) else: amap2D = amaps.UCLS(Mr, U) self.amap = np.reshape(amap2D, (h, w, U.shape[0])) self.m = h self.n = w self.q = U.shape[0] return self.amap def __str__(self): return 'pysptools.abundance_maps.amaps_int.UCLS object, maps: {0}x{1}x{2}'.format(self.m, self.n, self.q) @PlotInputValidation('UCLS') def plot(self, path, mask= None, interpolation='none', colorMap='jet', columns=None, suffix=None): _plot_abundance_map(self.amap, path, 'UCLS', mask, interpolation, colorMap, columns, suffix) @DisplayInputValidation('UCLS') def display(self, mask= None, interpolation='none', colorMap='jet', columns=None, suffix=None): _plot_abundance_map(self.amap, None, 'UCLS', mask, interpolation, colorMap, columns, suffix)
_document(UCLS)
[docs]class NNLS(object): """ NNLS performs non-negative constrained least squares with the abundance nonnegative constraint (ANC). Utilizes the method of Bro. """ def __init__(self): self.amap = None self.m = None self.n = None self.q = None @MapInputValidation('NNLS') def map(self, M, U, normalize=False, mask=None): """ NNLS performs non-negative constrained least squares of each pixel in M using the endmember signatures of U. Parameters: M: `numpy array` A HSI cube (m x n x p). U: `numpy array` A spectral library of endmembers (q x p). normalize: `boolean [default False]` If True, M and U are normalized before doing the spectra mapping. mask: `numpy array [default None]` A binary mask, when *True* the selected pixel is unmixed. Returns: `numpy array` An abundance maps (m x n x q). """ h,w,numBands = M.shape if normalize == True: M = _normalize(M) U = _normalize(U) Mr = np.reshape(M, (w*h, numBands)) if isinstance(mask, np.ndarray): m = np.reshape(mask, (w*h)) cMr = _compress(Mr, m) c_amap2D = amaps.NNLS(cMr, U) amap2D = _expand(c_amap2D, m, w*h, U.shape[0]) else: amap2D = amaps.NNLS(Mr, U) self.amap = np.reshape(amap2D, (h, w, U.shape[0])) self.m = h self.n = w self.q = U.shape[0] return self.amap def __str__(self): return 'pysptools.abundance_maps.amaps_int.NNLS object, maps: {0}x{1}x{2}'.format(self.m, self.n, self.q) @PlotInputValidation('NNLS') def plot(self, path, mask= None, interpolation='none', colorMap='jet', columns=None, suffix=None): _plot_abundance_map(self.amap, path, 'NNLS', mask, interpolation, colorMap, columns, suffix) @DisplayInputValidation('NNLS') def display(self, mask= None, interpolation='none', colorMap='jet', columns=None, suffix=None): _plot_abundance_map(self.amap, None, 'NNLS', mask, interpolation, colorMap, columns, suffix)
_document(NNLS)
[docs]class FCLS(object): """ Performs fully constrained least squares. Fully constrained least squares is least squares with the abundance sum-to-one constraint (ASC) and the abundance nonnegative constraint (ANC). """ def __init__(self): self.amap = None self.m = None self.n = None self.q = None @MapInputValidation('FCLS') def map(self, M, U, normalize=False, mask=None): """ Performs fully constrained least squares of each pixel in M using the endmember signatures of U. Parameters: M: `numpy array` A HSI cube (m x n x p). U: `numpy array` A spectral library of endmembers (q x p). normalize: `boolean [default False]` If True, M and U are normalized before doing the spectra mapping. mask: `numpy array [default None]` A binary mask, when *True* the selected pixel is unmixed. Returns: `numpy array` An abundance maps (m x n x q). """ h,w,numBands = M.shape if normalize == True: M = _normalize(M) U = _normalize(U) Mr = np.reshape(M, (w*h, numBands)) if isinstance(mask, np.ndarray): m = np.reshape(mask, (w*h)) cMr = _compress(Mr, m) c_amap2D = amaps.FCLS(cMr, U) amap2D = _expand(c_amap2D, m, w*h, U.shape[0]) else: amap2D = amaps.FCLS(Mr, U) self.amap = np.reshape(amap2D, (h, w, U.shape[0])) self.m = h self.n = w self.q = U.shape[0] return self.amap def __str__(self): return 'pysptools.abundance_maps.amaps_int.FCLS object, maps: {0}x{1}x{2}'.format(self.m, self.n, self.q) @PlotInputValidation('FCLS') def plot(self, path, mask= None, interpolation='none', colorMap='jet', columns=None, suffix=None): _plot_abundance_map(self.amap, path, 'FCLS', mask, interpolation, colorMap, columns, suffix) @DisplayInputValidation('FCLS') def display(self, mask= None, interpolation='none', colorMap='jet', columns=None, suffix=None): _plot_abundance_map(self.amap, None, 'FCLS', mask, interpolation, colorMap, columns, suffix)
_document(FCLS)