% note: This program runs on disguise experiment on AR database
% (partitioning part). For the occlusion experiment on Extended YaleB, you
% can just change some parameter's value.
%
% TrainDAT.mat  ----      the training data
%               ----      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 = 30  : the width of each block (AR)
%               ----      im_h = 21  : the height of each block (AR)
%                         The original image's size is 83*60, following SRC, we devide the image into 8 block
%                         (4*2), each with the size of about 21*30, and in our
%                         paper, we set the downsample gabor image with the
%                         size of 4*5 in each block.
%
% TestDAT.mat   ----      the testing data
%               ----      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;

im_h           =    21;            % for AR
im_w           =    30;            % for AR

par.ds_w       =    5;              % the downsample image's width in Gabor for AR
par.ds_h       =    4;              % the downsample image's heigth in Gabor for AR
par.ke_w       =    21;             % Gabor kernel's width
par.ke_h       =    21;             % 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_l   =    0.0005;         % parameter of l1_ls in learning
par.lambda_t   =    0.0005;         % parameter of l1_ls in testing
par.dim        =    20;            % the occlusion column dimension
par.dl_nIter   =    15;             % the maximal iteration number in occlusion dictionary learning

% this partition is just for the downsampled (83*60) image in AR database.
h_cor(1,:)     =   [1:21];       w_cor(1,:)     =   [1:30];
h_cor(2,:)     =   [22:42];      w_cor(2,:)     =   [1:30];
h_cor(3,:)     =   [43:63];      w_cor(3,:)     =   [1:30];
h_cor(4,:)     =   [63:83];      w_cor(4,:)     =   [1:30];
h_cor(5,:)     =   [1:21];       w_cor(5,:)     =   [31:60];
h_cor(6,:)     =   [22:42];      w_cor(6,:)     =   [31:60];
h_cor(7,:)     =   [43:63];      w_cor(7,:)     =   [31:60];
h_cor(8,:)     =   [63:83];      w_cor(8,:)     =   [31:60];

for  i = 1 : 8
ipts.dat    =    tr_dat(h_cor(i,:),w_cor(i,:),1:10);
[opts]      =    Create_GaborF(ipts, par); 
opts.gdat   =    opts.gdat./ repmat(sqrt(sum(opts.gdat.*opts.gdat)),[size(opts.gdat,1) 1]); % unit norm 2
tr(i).dat   =    opts.gdat;
end

for  i = 1 : 8
ipts.dat    =    tt_dat(h_cor(i,:),w_cor(i,:),1:10);
[opts]      =    Create_GaborF(ipts, par);
opts.gdat    =    opts.gdat./ repmat(sqrt(sum(opts.gdat.*opts.gdat)),[size(opts.gdat,1) 1]); % unit norm 2
tt(i).dat   =    opts.gdat;
end

occlu_dat      =    eye(im_h*im_w);
oclu_ipts.dat  =    uint8(reshape(occlu_dat,[im_h im_w im_h*im_w]));
[oclu_opts]    =    Create_GaborF (oclu_ipts, par);

D_input        =    randn(size(oclu_opts.gdat,1),par.dim); % Occlusion dictionary initialization
D_input        =    D_input./ repmat(sqrt(sum(D_input.*D_input)),[size(D_input,1) 1]); % unit norm 2
[D_Oclu error] =    Occlusion_Dict_Learn (oclu_opts.gdat,D_input,  par.dl_nIter, par.lambda_l);


classids       =    unique(tr_lab);
NumClass       =    length(classids);
tt_num         =    size(tt_dat,3);
tt_ID          =    [];

for js  =   1:8
    D   =  [tr(js).dat D_Oclu];
    
    for i  =  1: tt_num
    [s,status]   =  l1_ls(D, tt(js).dat(:,i), par.lambda_t); % l1_ls toolbox
    gap          =  zeros(1,NumClass);
    
        for j   =  1:NumClass
        temp_s      =  zeros(size(tr(js).dat,2),1);
        temp_s(find(j==tr_lab))=s(find(j==tr_lab));
        add_coef    =   s(size(tr(js).dat,2)+1:end,1);
        zz          =   tt(js).dat(:,i)-tr(js).dat *temp_s-D_Oclu*add_coef;
        gap(j)      =   zz(:)'*zz(:); 
        end
    
    index           =  find(gap==min(gap));
    tt_ID(js,i)     =  classids(index(1));
    end
    
end

indentify     =    [];

for i  =  1 :size(tt_num,2)
       temp    =    tt_ID(:,i);
       
       for j_s=1:8
           N(j_s)  =  length(find(temp==temp(j_s)));
       end
       
       index         =   find(N==max(N));
       indentify(i)  =  temp(index(1));
end

reco_ratio     =     (sum(indentify==tt_lab))/length(tt_lab);