% This main function is related to Face Recognition without occlusion by
% using GSRC of our paper in ECCV 2010.
% note:
% TrainDAT.mat  ----      the training data of Extended YaleB (or AR or FERET pose subset)
%               ----      tr_dat: 3d Training data matrix, each tr_dat(:,:,i) is an (downsampled) 
%                            image matrix.
%               ----      tr_lab: a row vector, each element is the label of each training sample
%               ----      im_w  : the width of image
%               ----      im_h  : the height of image
%
% TestDAT.mat   ----      the testing data of Extended YaleB (or AR or FERET pose subset)
%               ----      tt_dat: 3d Training data matrix, each tr_dat(:,:,i) is an (downsampled) 
%                            image matrix.
%               ----      tt_lab: a row vector, each element is the label
%               of each training sample
%
% Copyright. Mike Yang, PolyU, Hong Kong.
% reference:
% M. Yang and L. Zhang.  Gabor feature based sparse representation
% for face recognition with Gabor occlusion dictionary. In ECCV 2010.
% S.J. Kim, K. Koh, M. Lustig, S. Boyd, and D. Gorinevsky. A method
% for large-scale l1-regularized least squares. IEEE Journal on Selected Topics in Signal Processing, 1(4):606¨C617, 2007.

clear all
clc

path      =     cd;
addpath([path '\utilities']);

load TrainDAT.mat;
load TestDAT.mat;

tr_ipts.dat    =    tr_dat;
tt_ipts.dat    =    tt_dat;

par.ds_w       =    floor(im_w/8.1);      % the downsample image's width in Gabor
par.ds_h       =    floor(im_h/8);        % the downsample image's heigth in Gabor
par.ke_w       =    31;                   % Gabor kernel's width
par.ke_h       =    31;                   % Gabor kernel's heigth
par.raT        =    0.9;                  % Gabor kernel's energy preseving ratio
par.Kmax       =    pi/2;                 % Gabor kernel's para, default(pi/2)
par.f          =    sqrt(2);              % Gabor kernel's para, default(sqrt(2))
par.sigma      =    pi;                   % Gabor kernel's para, default(pi or 1.5pi)
par.Gabor_num  =    40;                   % Gabor kernel's number, default (5 scales and 8 orientations)
par.lambda     =    0.001;                % parameter of l1_ls in testing
par.dim        =    30;                   % the pca dimension

% create Gabor feature if you don't have
[tr_opts]      =    Create_GaborF (tr_ipts, par);
train_temp     =    tr_opts.gdat;
save('GaborTrain','train_temp'); clear tr_opts
[tt_opts]      =    Create_GaborF (tt_ipts, par);
test_temp      =    tt_opts.gdat;
save('GaborTest','test_temp');   clear tt_opts

% load Gabor feature if you have had
% load GaborTrain;
% load GaborTest;

[disc_set,disc_value,Mean_Image]  =  Eigenface_f(train_temp,par.dim);  

tr_gdat        =    disc_set'* train_temp;
tt_gdat        =    disc_set'* test_temp;
tr_gdat        =    tr_gdat./ repmat(sqrt(sum(tr_gdat.*tr_gdat)),[size(tr_gdat,1) 1]); % unit norm 2
tt_gdat        =    tt_gdat./ repmat(sqrt(sum(tt_gdat.*tt_gdat)),[size(tt_gdat,1) 1]); % unit norm 2

classids       =    unique(tr_lab);
NumClass       =    length(classids);
tt_num         =    size(tt_dat,3);
tt_ID          =    [];
    
for i  =  1: tt_num
    
    [s,status] = l1_ls(tr_gdat, tt_gdat(:,i), par.lambda); % l1_ls toolbox
    gap        = zeros(1,NumClass);
    
    for j   =  1:NumClass
        temp_s =  zeros(size(s));
        class  =  classids(j);
        temp_s(find(class==tr_lab))  =  s(find(class==tr_lab));
        zz     =  tt_gdat(:,i)-tr_gdat*temp_s;
        gap(j) =  zz(:)'*zz(:); 
    end
    
    index  =  find(gap==min(gap));
    tt_ID  =  [tt_ID classids(index(1))];
end

reco_ratio=(sum(tt_ID==tt_lab))/length(tt_lab); % output the recognition rate
