www.pudn.com > gene-selection.zip > TRIP.m, change:2012-08-10,size:2232b


%10th August 2012 
%===================================================================== 
% Programmer:  
% ArulMurugan Ambikapathi, E-mail: a.arulmurugan@gmail.com 
% ------------------------------------------------------- 
% Reference:  
% A.Ambikapathi, T.-H. Chan, C.-Y. Chi, and  K. Keizer, ``Hyperspectral 
% data geometry based estimation of number of endmembers using p-norm 
% based pure pixel identification algorithm,'' to appear 
% in IEEE Trans. Geoscience and Remote Sensing. 
%====================================================================== 
% An implementation of p-norm purest pixel identification (TRI-P) Algorithm 
% [PP_Index]=TRIP(X,N,p,N_max,CNN) 
%====================================================================== 
%  Input 
%  X is M-by-L data matrix where M is the spectral bands (or observations) and L is the number of pixels (data length).    
%  N is the current iteration number in GENE-CH or GENE-AH 
%  p denotes the pth norm say, p=2 
%  N_max is maximum bound on the number of endmembers, where N <= N_max <= M, say N_max=25. 
%  CNN is an estimate of the noise covariance matrix. 
%---------------------------------------------------------------------- 
%  Output 
%  PP_Index is the set of N purest pixel indices 
%======================================================================== 
 
function [PP_Index]=TRIP(X,N,p,N_max,CNN)      
Rw=CNN; 
[M,L] = size(X);  
d = mean(X,2); 
U = X-d*ones(1,L); 
R = U*U'; 
Rx_true = (R./L) - Rw; 
[eV D] = eig(Rx_true); 
C = eV(:,end:-1:M-(N_max-2)); 
Xd_t = C'*(X-d*ones(1,L));  % Dimension reduced data vectors (Xd is (N-1)-by-L) 
 
%===algorithm========= 
A_set=[]; Xd = [Xd_t; ones(1,L)]; index = [];               % N is endmember, p is p-norm, Xd is N by L 
[val ind]=max(sum(abs(Xd).^p).^(1/p));    
A_set = [A_set Xd(:,ind)]; 
index = [index ind];                                        % save index 
for i=2:N     
    XX = (eye(N_max) - A_set * pinv(A_set)) * Xd;     
    [val ind]=max(sum(abs(XX).^p).^(1/p));                  %    N*p'(p-norm)*N* 
    A_set = [A_set Xd(:,ind)];                              % find A_est one by one 
    index = [index ind];                                    % save index 
end 
PP_Index=index;