Source code for pysptools.material_count.vd

# 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
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# See the License for the specific language governing permissions and
# limitations under the License.
# - This file is part of the PySptools package.

HfcVd function

from __future__ import division

import numpy as np
import scipy as sp
import scipy.stats as ss

def est_noise(y, noise_type='additive'):
    This function infers the noise in a
    hyperspectral data set, by assuming that the
    reflectance at a given band is well modelled
    by a linear regression on the remaining bands.

        y: `numpy array`
            a HSI cube ((m*n) x p)

       noise_type: `string [optional 'additive'|'poisson']`

    Returns: `tuple numpy array, numpy array`
        * the noise estimates for every pixel (N x p)
        * the noise correlation matrix estimates (p x p)

        Jose Nascimento ( and Jose Bioucas-Dias (
        For any comments contact the authors
    def est_additive_noise(r):
        small = 1e-6
        L, N = r.shape
        w=np.zeros((L,N), dtype=np.float),r.T)
        RRi = np.linalg.pinv(RR+small*np.eye(L))
        RRi = np.matrix(RRi)
        for i in range(L):
            XX = RRi - (RRi[:,i]*RRi[i,:]) / RRi[i,i]
            RRa = RR[:,i]
            RRa[i] = 0
            beta =, RRa)
            w[i,:] = r[i,:] -,r)
        Rw = np.diag(np.diag(,w.T) / N))
        return w, Rw

    y = y.T
    L, N = y.shape
    #verb = 'poisson'
    if noise_type == 'poisson':
        sqy = np.sqrt(y * (y > 0))
        u, Ru = est_additive_noise(sqy)
        x = (sqy - u)**2
        w = np.sqrt(x)*u*2
        Rw =,w.T) / N
    # additive
        w, Rw = est_additive_noise(y)
    return w.T, Rw.T

[docs]def hysime(y, n, Rn): """ Hyperspectral signal subspace estimation Parameters: y: `numpy array` hyperspectral data set (each row is a pixel) with ((m*n) x p), where p is the number of bands and (m*n) the number of pixels. n: `numpy array` ((m*n) x p) matrix with the noise in each pixel. Rn: `numpy array` noise correlation matrix (p x p) Returns: `tuple integer, numpy array` * kf signal subspace dimension * Ek matrix which columns are the eigenvectors that span the signal subspace. Copyright: Jose Nascimento ( & Jose Bioucas-Dias ( For any comments contact the authors """ y=y.T n=n.T Rn=Rn.T L, N = y.shape Ln, Nn = n.shape d1, d2 = Rn.shape x = y - n; Ry =, y.T) / N Rx =, x.T) / N E, dx, V = np.linalg.svd(Rx) Rn = Rn+np.sum(np.diag(Rx))/L/10**5 * np.eye(L) Py = np.diag(,,E))) Pn = np.diag(,,E))) cost_F = -Py + 2 * Pn kf = np.sum(cost_F < 0) ind_asc = np.argsort(cost_F) Ek = E[:, ind_asc[0:kf]] return kf, Ek # Ek.T ?
# Comments on using complex number: # # Use only scipy and numpy functions for a correct use of complex number. # # scipy.sqrt() deal by the book with complex number, # it's more tricky when using math and numpy modules. #
[docs]def HfcVd(M, far='default'): """ Computes the vitual dimensionality (VD) measure for an HSI image for specified false alarm rates. When no false alarm rate(s) is specificied, the following vector is used: 1e-3, 1e-4, 1e-5. This metric is used to estimate the number of materials in an HSI scene. Parameters: M: `numpy array` HSI data as a 2D matrix (N x p). far: `list [default default]` False alarm rate(s). Returns: python list VD measure, number of materials estimate. References: C.-I. Chang and Q. Du, "Estimation of number of spectrally distinct signal sources in hyperspectral imagery," IEEE Transactions on Geoscience and Remote Sensing, vol. 43, no. 3, mar 2004. J. Wang and C.-I. Chang, "Applications of independent component analysis in endmember extraction and abundance quantification for hyperspectral imagery," IEEE Transactions on Geoscience and Remote Sensing, vol. 44, no. 9, pp. 2601-1616, sep 2006. """ N, numBands = M.shape # calculate eigenvalues of covariance and correlation between bands lambda_cov = np.linalg.eig(np.cov(M.T))[0] # octave: cov(M') lambda_corr = np.linalg.eig(np.corrcoef(M.T))[0] # octave: corrcoef(M') # not realy needed: lambda_cov = np.sort(lambda_cov)[::-1] lambda_corr = np.sort(lambda_corr)[::-1] if far == 'default': far = [10**-3, 10**-4, 10**-5] else: far = [far] numEndmembers_list = [] for y in range(len(far)): numEndmembers = 0 pf = far[y] for x in range(numBands): sigmaSquared = (2.*lambda_cov[x]/N) + (2.*lambda_corr[x]/N) + (2./N)*lambda_cov[x]*lambda_corr[x] sigma = sp.sqrt(sigmaSquared) tau = -ss.norm.ppf(pf, 0, abs(sigma)) if (lambda_corr[x]-lambda_cov[x]) > tau: numEndmembers += 1 numEndmembers_list.append(numEndmembers) return numEndmembers_list