function [U] = hyperNfindr(M, q)
% HYPERNFINDR Performs the N-FINDR (endmember extraction) algorithm
% Performs the N-FINDR algorithm to find q endmembers. If only M is
% given as input, this function calls hyperHfcVd to estimate the number
% of endmembers (q) and then hyperPct to reduce dimensionality to (q-1).
%
% Usage
% [U] = hyperNfindr(M)
% [U] = hyperNfindr(M, q)
% Inputs
% M - 2d matrix of HSI data (p x N)
% q - Number of endmembers to find
% -- if not given, q is obtained from hyperHfcVd(M, 10^-3)
% Outputs
% U - Recovered endmembers (p x q)
%
% References
% M. Winter, "N-findr: an algorithm for fast autonomous
% spectral endmember determination in hyperspectral data," SPIE’s
% International Symposium on Optical Science, Engineering, and
% Instrumentation, pages 266–275. International Society for Optics
% and Photonics, 1999.
% Error trapping
if ndims(M) ~= 2
warning('WarnTests:dim', ...
'Input image must be p x N.\n',...
'Converting with hyperConvert2d.\n')
M = hyperConvert2d(M);
end
M_orig = M;
[p, N] = size(M);
if nargin == 1
fprintf('Implementing hyperHfcVd to determine the number of endmembers.\n')
q = hyperHfcVd(M_orig, [10^-3]);
fprintf('Reducing dimensionality to (q-1) using hyperPct.\n')
M = hyperPct(M, q-1);
elseif q < p+1
warning('WarnTests:dim', ...
strcat('N-FINDR requires (q-1) spectral bands.\n',...
'Performing PCA to reduce dimensionality.\n'))
M = hyperPct(M, q-1);
elseif q > p+1
warning('WarnTests:dim', ...
strcat('N-FINDR requires (q-1) spectral bands.\n',...
'Performing PCA to reduce dimensionality.\n'))
error('ErrTests:dim', ...
strcat('N-FINDR cannot find more than (p+1) endmembers (q),\n', ...
'where p is the number of available spectral bands.\n'))
end
% Initialize
U_idx = randperm(N,q); % Random endmember selection
E = M(:,U_idx); % Endmember matrix
V = abs(det([ones(1,q); E])) / factorial(q-1); % Simplex volume
vols = zeros(q,1);
% Search for maximum volume simplex
for j = 1:N;
% Replace each column of E with sample vector M(:,j)
% and compute the volume for each
for k = 1:q;
E_tmp = E;
E_tmp(:,k) = M(:,j);
vols(k) = abs(det([ones(1,q); E_tmp])) / factorial(q-1);
end
% If max volume is greater than previous V, update E and V
[V_tmp,k_idx] = max(vols);
if V_tmp > V
V = V_tmp;
E(:,k_idx) = M(:,j);
U_idx(k_idx) = j;
end
end
% Return endmembers
U = M_orig(:, U_idx);
评论0