function [axes,Pcomponents] = twodlpp2(images,label,NN,numP)
%%  function [axes,Pcomponents] = twodlpp2(images,label,NN,numP)
%%     It is to compute the first several principal axes 
%%  and the correspondent principal components of 2-dimensional LPP.
%%  Input:   images:   (M x N x K) matrix containing K images
%%                     which sizes are all (M x N).
%%           label :   (1 x K) numeric vector (i.e. 1,2,...,c) indicating 
%%                     which class each point is belong to. 
%%           NN    :   number of neighbor
%%           numP  :   the number of principal axes specified.
%%  Output:  axes  :   (N x numP) matrix, the i'th colomn of which
%%                     is the i'th principal axis.
%%           Pcomponents: (M x numP x K) matrix, containing K feature
%%                         matrix with size (M x numP),whose i'th colomn 
%%                         is the i'th principal component.
%%
%%  coded by sbchen, March 2005.

[M,N,K] = size(images);
%images = double(images);

S = nn_graph_c(reshape(images,M*N,K),label,NN,1);
%S = nn_graph_c(reshape(images,M*N,K),label,NN,0);

D = diag(sum(S));
L = D-S;              %Laplacian matrix

% % DD = kron(D,eye(M));   %out of memory in 1G memory
% % LL = kron(L,eye(M));
% % A=zeros(0,N);
% % for i=1:K
% %     A = [A;images(:,:,i)];
% % end
% % [axes,dd]=eig(A'*LL*A,A'*DD*A);
prod = images(:,:,K)'*images(:,:,K);
DD = D(K,K)*prod;
LL = L(K,K)*prod;
for i=1:(K-1)
    prod = images(:,:,i)'*images(:,:,i);
    DD = DD + D(i,i)*prod;
    LL = LL + L(i,i)*prod;
    for j=(i+1):K
        prod = L(i,j)*images(:,:,i)'*images(:,:,j);
        LL = LL + prod + prod';
    end
end

[axes,dd]=eig(LL,DD);

% cD = diag(dd);
% [cD,I] = sort(cD);
% idx = I(find(cD > 0));
% axes = U(:,idx);

if nargin==4 
    axes = axes(:,1:numP);
end

for i = 1:K
    Pcomponents(:,:,i) = images(:,:,i)*axes;
end