function [im2, d]  = mv_pca(im1, neig)

% function imeig = pca(im1, neig)
%
% principle component analysis of image block
% assumes appropriate means already removed (see stripmean)
%
% im1   = input ni x nj x nt array
% neig  = number of modes to keep
%
% im2  = output ni x nj x neig array

[ni, nj, nt] = size(im1);
im1 = reshape(im1, ni*nj, nt);

if ni*nj < nt
   S = im1*im1';
else
   S = im1'*im1;
end

[evec, eval] = eig(S);
d=diag(eval);
f(1); clf; semilogy(d);
[d, i] = sort(-d);
d = -d;
evec = evec(:, i);
evec = evec(:, 1:neig);

if ni*nj < nt
   im2 = reshape(evec, ni, nj, neig);
else
   im2 = im1*evec;
   for i = 1:neig
      s = sqrt(im2(:, i)'*im2(:, i));
      im2(:, i) = im2(:, i)/s;
   end
   im2 = reshape(im2, ni, nj, neig);
end