Endmember Extraction Using the Pixel Purity Index (PPI) Algorithm

Endmember Extraction Using the Pixel Purity Index (PPI) Algorithm The Pixel Purity Index (PPI) algorithm is used to extract endmember signatures from hyperspectral data cubes. The algorithm identifies endmember pixels by projecting dimensionality-reduced data onto randomly generated unit vectors and finding the pixels with the maximum projection values.

Algorithm Principle

Perform dimensionality reduction on the hyperspectral data cube using Principal Component Analysis (PCA), Maximum Noise Fraction (MNF), or no reduction at all. Generate a set of random unit vectors (default 10,000). Project the reduced data onto each random unit vector and record the maximum projection value for each pixel. Rank the pixels based on their maximum projection values and select the top N pixels as the endmember signatures. ## Algorithm Implementation Here's how to implement the PPI algorithm in MATLAB:

   function endmemSig = ppi(cube, num, option)
%PPI Extract endmember signatures using pixel purity index.
%   ENDMEMSIG = PPI(CUBE, NUM) estimates the endmember signatures by
%   projecting a reduced data on randomly generated unit vectors. NUM
%   describes the number of endmember signatures present in the CUBE. CUBE
%   is a numeric array with dimensions M-by-N-by-C containing hyperspectral
%   data or a hypercube object. M and N are rows and columns of the
%   hyperspectral cube respectively. C describes number of bands in the
%   hyperspectral cube. ENDMEMSIG is a numeric matrix and it provides the
%   endmember signatures of dimensions C-by-NUM.
%
%   ENDMEMSIG = PPI(___, Name, Value) estimates the endmember signatures by
%   projecting a randomly generated unit vectors on a reduced data with
%   selected NV pair.
%
%   'NumVectors'    Number of randomly generated unit vectors. As the
%                   number increases algorithm does a better job of finding
%                   the endmembers along with the increasing the
%                   computation complexity. Default value is 10^4.
%
%   'ReductionMethod': Method used to perform dimension reduction. Choices
%   are:
%
%       'PCA'  Based on the maximum data variation Principal Component
%              Analysis (PCA) transforms large set of data into smaller one
%              by maintaining most of the information.
%
%       'MNF'  Based on the Noise fractions Maximum Noise Fraction (MNF)
%              transforms large set of data into smaller one by maintaining
%              most of the information. MNF is a method which orders the
%              components according to image quality. The default method is
%              MNF.
%
%       'NONE' Project original data on randomly generated unit vectors.
%
%   Notes
%   -----
%   1) PPI algorithm uses default random seeds for generating random unit
%   vectors.
%
%   Class Support
%   -------------
%   The input array CUBE must be one of the following classes: uint8,
%   uint16, uint32, int8, int16, int32, single, double, uint64, int64 or
%   hypercube object. It must be real and non-sparse. NUM and NumVectors
%   are positive integer-valued numeric scalar. The output ENDMEMSIG is a
%   numeric matrix of same class as input cube with size C-by-NUM.
%
%   Reference
%   ---------
%   J.W  Boardman, F.A. Kruse and  R.O. Green, "Mapping target signatures
%   via partial unmixing of AVIRIS data", Technical Report, California,
%   USA, 1995.
%
%   Example - 1 : Find endmember signatures.
%   -----------
%  % Load indian pines dataset
%   cube = load('indian_pines.mat');
%
%   % Estimate number of endmembers.
%   num = countEndmembersHFC(cube.indian_pines);
%
%   % Estimate and visualize the endmember signatures
%   endmemSig = ppi(cube.indian_pines, num);
%   plot(endmemSig);
%   xlabel('Band index')
%   ylabel('Signatures')
%   title('Endmembers')
%
%   Example - 2 : Find endmember signatures with 'pca' dimension reduction.
%   % method.
%   -----------
%   % Load hyperspectral cube.
%   obj = hypercube('paviaU.hdr');
%
%   % Estimate number of endmembers.
%   num = countEndmembersHFC(obj);
%
%   % Estimate and visualize the endmember signatures using PPI.
%   endmemSig = ppi(obj, num, 'ReductionMethod', 'pca', 'NumVectors', 10^3);
%   plot(endmemSig);
%   xlabel('Band index')
%   ylabel('Signatures')
%   title('Endmembers')
%
%   See also nfindr, fippi, hypercube, countEndmembersHFC.

%  Copyright 2020-2021 The MathWorks, Inc.


arguments
    cube {mustBeA(cube,'hypercube')}
    num (1,1) {mustBeNumeric, mustBeNonNan, mustBeFinite, mustBeNonsparse, ...
        mustBeNonempty, mustBeReal, mustBeInteger, mustBePositive, ...
        mustBeLessSize(num, cube)}
    option.NumVectors (1,1) {mustBeNumeric, mustBeNonNan, mustBeFinite, ...
        mustBeNonsparse, mustBeNonempty, mustBeReal, mustBeInteger,...
        mustBePositive} = 10^4
    option.ReductionMethod (1,1) string ...
        {validatestring(option.ReductionMethod,{'MNF', 'PCA', 'NONE'})} = "MNF"

end

% For compiler tests
if isdeployed
    rootDir = ctfroot;
else
    rootDir = matlabroot;
end
% Register resources/hyperspectral.
matlab.internal.msgcat.setAdditionalResourceLocation(rootDir);

% Input parsing.
if isobject(cube)
    cube = cube.DataCube;
end

numOfUnitVec = option.NumVectors;
method = validatestring(option.ReductionMethod, {'MNF', 'PCA', 'NONE'});

% Volume dimensions.
[rows, cols, channels] = size(cube);
numOfSamples = rows*cols;

originalCube = reshape(cube, numOfSamples, channels);

% Cast the input data to single data type
if isinteger(cube)
    cube = single(cube);
    num = single(num);
end

if strcmpi(method, 'PCA')
    % Principal Component Analysis (PCA) on hyperspectral cube.
    cube = reshape(hyperpca(cube, num),numOfSamples, num);

elseif strcmpi(method, 'MNF')
    % Maximum noise fraction (MNF) transformation is used to reduce the
    % computational requirement.
    cube = reshape(hypermnf(cube, num), numOfSamples, num);

else
    cube = reshape(cube, numOfSamples, channels);
end

cube = cube';

% Determine the maximum projection.
maxProjections = cast(zeros(size(cube, 2), 1), class(cube));

% Generate random unit vectors.
prevS = rng('default');
unitVec = randn(numOfUnitVec, size(cube,1));

% Clean up the random state.
cleanup = onCleanup(@() rng(prevS));
for numVect = 1:numOfUnitVec
    [~, idx] =  max(abs(unitVec(numVect,:) *cube));
    maxProjections(idx) = maxProjections(idx) + 1;
end

% Find the maximum projected votes to decide the endmember signatures.
[~, idx] = sort(maxProjections, 'descend');

% Extract endmembers form the volume.
endmemSig = originalCube(idx(1:num),:)';
end

function mustBeA(hypercube,classIn)

if isobject(hypercube)
    validateattributes(hypercube,{classIn},{'nonempty'}, 'ppi','hypercube');
    cube = hypercube.DataCube;
else
    cube = hypercube;
end
validateattributes(cube, ...
    {'numeric'}, {'nonsparse', 'nonempty', 'real','nonnan', ...
    'finite', 'ndims', 3}, mfilename, 'hypercube', 1);

end

% Validate the numComponents are less than number of bands.
function mustBeLessSize(a, b)
if isobject(b)
    cube = b.DataCube;
else
    cube = b;
end
sz = size(cube);
validateattributes(a, {'numeric'},...
    {'>', 0, '<=',sz(3),'scalar'},mfilename, 'num',2);

if (a >= sz(1)*sz(2))
    error(message('hyperspectral:hyperpca:notGreaterEqual'));
end
end

In summary, the PPI algorithm is a geometry-based unsupervised endmember extraction method that identifies endmember pixels by projecting dimensionality-reduced data onto random unit vectors and finding the pixels with the maximum projection values. It is commonly used for endmember extraction and nonlinear unmixing analysis of hyperspectral images.