% Principal Component Analysis 
% y=W.'(x-m)
% INPUT
% X : d*n matrix. (d>>n)
%   X(:,i) is a feature vector of d dimension 
%   n : number of samples 
% k : number of principal components to be produced.
% OUTPUT
% model.
%       W : d*k matrix, containing k largest eigenvectors of cov(x)
%       m : mean of X.
% NOTE
% 1.Input X can be obtained using the following code:
% for i=1:n
%     x=imread(filename); 
%     X(:,i)=x(:);
% end
% 2. To display the eigenvector as an image (like eigenface), please refer
% to the following code:
%  im=reshape(W(:,i),sz(1),sz(2)); //sz the original size of image.
%  subplot(...);    
%  imagesc(im);
%  colormap(gray(255)); //adjust the colormap to gray(255) i.e. 0~255 gray value
%  title(sprintf('eignface %d',i));
% 3. To project a new feature vector to PCA space, refer to the following
% code
%   y=model.W.' *(x-model.m);    //x is the feature vector
%
% 4. Reconstruct from PCA space;
%   x1=model.W * y + model.m;   
%
% 5. Be careful of variable type casting.
%      This function used double(.) to cast integers to double
%      If you prefer im2double to normalize image pixel values to 0..1,
%      please change accordingly.
%
% by guodong (http://guod.cn)
% last update: 2008-2-26

function model=pca(X,k)

[d n]=size(X);
X=im2double(X);
m=mean(X,2);
X=X-m*ones(1,n);

[V,D]=eig(X.'*X);
V=fliplr(V);
D=flipud(diag(D));
V=X*V*diag((D).^(-1/2));  % normalize V

model.W=V(:,1:k);
model.m=m;
model.D=D;
model.k=k;

return;