%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Multi-Variate Gaussian Classifier
% demo for COMP 435
% Z. Li, 2010.02.10
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function []=face_classify()
doEigenFace = 1;
doFisherFace = 1;
doPlot = 1;
model = 'pca';  % 'pca' or 'lda'
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% (1) Eigenface - PCA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if (doEigenFace)
% download the att_face.zip to the directory below
facepath='D:/zli/teaching/comp435/matlab/data/att_faces/';
nSubj=40; nPic=10;
% icon size
h=24; w=20;
%count 
if (0)
t=1;
% create face icons 
% load face images, re-size to icons of size h x w
    for j=1:nSubj
        for k=1:nPic
            fim_name=sprintf('%ss%d/%d.pgm', facepath, j, k);
            im = imread(fim_name);
            icn = imresize(im, [h, w], 'bilinear');
            faces(t,1:w*h) = (icn(:))'; t=t+1;
        end
    end
    if (1)
        %save face_icons_20x18.mat faces;
        save face_icons_24x20.mat faces;
    end
else % 2nd time load it directly
    %load face_icons_20x18.mat;
    load face_icons_24x20.mat;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% pca of faces
% output
%   EigFaces - PCA basis
%   EigVal   - PCA eig values
%   indx     - sorted by eig values
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load face_icons_24x20.mat;

faces = double(faces);
mfaces = mean(faces); mfaces2d = reshape(mfaces, h, w);
% covariance
S = cov(faces);
% eigen value decomposition 
[EigFaces, L]=eig(S);
% notice that pca eigen values are sorted from small to big by default,
% have to correct that. 
[EigVal, indx] = sort(diag(L), 'descend');

figure(1); stem(EigVal(1:180), '.'); title('PCA eigen values');
figure(2);
subplot(2,4, 1); colormap('gray'); imagesc(mfaces2d); title('mean f_0'); 
% pca dimensions we want 
m = 7;
for k=1:m
    subplot(2, 4, k+1);
    eigf = EigFaces(:, indx(k));
    eigf2 = reshape(eigf, h, w); 
    if (1)
        colormap('gray'); imagesc(eigf2); axis off;
    else
        stem(eigf, '.');
    end
    str=sprintf('eigen f_%d', k'); title(str);
end

end % if (doEigenface)
% PCA - A1
kDim = 40; 
A1=double(EigFaces(:, indx([1:kDim])));


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% (2) Fisherface - LDA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% prepare labels: faceIds
nSubj=40; nImages = 10;
faceIds=zeros(nSubj*nImages, 1);
for k=1:nSubj
    offs=(k-1)*nImages;
    faceIds(offs+1:offs+nImages) = k;
end

if (doFisherFace)
    [A2, nClass]=myLDA(faces*A1, faceIds);
    m = 8; A3 = A1*A2; 
    figure(3); 
    for k=1:m
        subplot(2,4,k);
        ff = A3(:,k);
        ff2 = reshape(ff, h, w); 
        if (1)
            colormap('gray'); imagesc(ff2); axis off;
        else
            stem(ff, '.');
        end
        str=sprintf('fisher f_%d', k'); title(str);
    end
end %if (doFisherFace)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% (3) Projection
%  kDim - how many dimension to use
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% faces projected to x:
if model == 'pca'
    A = A1(:,1:kDim);
else
    A = A1*A2;
end
% faces of R^(hw) are now projected to x of R^kDim
x = double(faces)*A;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% (3) Visualization
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% visualization 
plotOption = '2d';  % or '2d' 
if (doPlot)
    figure(4); grid on; hold on;  ylabel('x_2'); xlabel('x_1');
    str = sprintf('projected %s faces', model);  title(str);
    subj = [5 12 20 32 ]; pics =[1: 8]; 
    styls{1}='+r'; styls{2}='*g'; styls{3}='ob'; styls{4}='xm'; styls{5}='^y';
    % plot all faces
    if (0)
        if plotOption == '3d'
            plot3(x(:,1), x(:,2), x(:,3), '.k'); hold on;
        else   % '2d'
            plot(x(:,1), x(:,2), '.k'); hold on;
        end
    end
    
    for k=1:length(subj)
        offs = 10*(subj(k)-1);
        if plotOption == '3d'
            plot3(x(offs+1:offs+5, 1), x(offs+1:offs+5, 2),x(offs+1:offs+5, 3), styls{k});
        else
            plot(x(offs+1:offs+5, 1), x(offs+1:offs+5, 2), styls{k});
        end
    end
end % plot

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% (3) Classifiers
%   kNN 
%   Bayesian
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% query face ids
qSubj=[5 12 20 32 38]; 
% query face image number
qImg = 7;
% pull query faces: qx
for k=1:length(qSubj)
        qOffs(k) = 10*(qSubj(k)-1)+qImg;
        qx(k,:)=x(qOffs(k), :);
end
% training samples: tx
% setdiff(a, b), a must contain b !
% find non-overlapping training set:
txOffs = setdiff([1:nSubj*nImages], qOffs);
% training data and label in tx and tx_lables
tx = x(txOffs,:);
tx_labels = faceIds(txOffs); 

%%%%%%%%%%%%%%%%%%%%%%
% test classifiers
%%%%%%%%%%%%%%%%%%%%%%
qDim = 3; qx = qx(:,1:qDim); tx = tx(:,1:qDim); 
[rec_labels1, err]=classify(qx, tx, tx_labels, 'linear');
[rec_labels2, err]=classify(qx, tx, tx_labels, 'quadratic');
[rec_labels3] = knnClassify(qx, tx, tx_labels); 
fprintf('\n Probe face ids: '); fprintf('%d ', qSubj); 
fprintf('\n Recognized ids, Bayesian-linear : '); fprintf('%d ', rec_labels1); 
fprintf('\n Recognized ids, Bayesian-quadratic : '); fprintf('%d ', rec_labels2); 
fprintf('\n Recognized ids, 1-NN classify : '); fprintf('%d ', rec_labels3); 


return;