function pca_plot(D, classes)

%
% pca_plot(D, [classes])
%
%  D          data matrix (n x dim)
%  [classes]  cluster id for each data vector (optional)
%

% 
% (C) Esa Alhoniemi, esa.alhoniemi@hut.fi
%

[l dim] = size(D);
comps   = 1:dim;

n_comps = length(comps);

% compute principal components

D_zeromean = D - ones(l,1) * mean(D);  
[V, S]     = eig(D_zeromean' * D_zeromean/(l-1));
 
% sort them in descending order
 
[y,ind] = sort(diag(S)); ind = flipud(ind);
V       = V(:,ind);

% process cluster id numbers

if nargin==2
  u = unique(classes);
  lu = length(u); 
else
  lu = 1;
  u  = 1;
  classes = ones(l,1);
end

clf

% set colormap

co = [[0 0 0];
      [1 0 0];
      [0 1 0];
      [0 0 1];
      [1 1 0];
      [1 0 1]];
	  	  
colormap(co);

% draw

n = 0;
for i = comps
  for j = comps
    n = n + 1;
    subplot(n_comps, n_comps, n); hold on
    if j > i
      for k = 1:lu
	f = find(classes == u(k));
	plot(D(f,:) * V(:,i), D(f,:) * V(:,j), '.', 'Markersize', ...
	     5, 'color', co(k,:));
      end
    elseif j==i
      axis([0 2 0 2])
      text(1, 1, sprintf('PC %d', i), 'horizontalalignment','center',...
	   'verticalalignment','middle'); 
    end
    axis off
  end
end