orl_src.asv

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% SRC, Sparse Representation-based Classification, 
% an algorithm in Pattern Recognition area
% 
% Author:
%   Written by Denglong Pan,
%   pandenglong@gmail.com
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc
t0=clock;
Image_row_NUM=112;Image_column_NUM=92; 
NN=Image_row_NUM*Image_column_NUM;%

Class_Train_NUM=5;

% total
Class_Sample_NUM=10; 
Class_Test_NUM=Class_Sample_NUM-Class_Train_NUM;

Class_NUM=40;
Train_NUM=Class_NUM*Class_Train_NUM; % 
Test_NUM=Class_NUM*(Class_Sample_NUM-Class_Train_NUM); % 

% Number of PCA principal component
Eigen_NUM=80;
Disc_NUM=Eigen_NUM;

% for Kernel Fisherface and CKFD regular and irregular
Dim_Begin=10; Dim_End=Disc_NUM; Dim_Interval=1; 

% is same as (Dim_End-Dim_Begin+Dim_Interval)/Dim_Interval
Dim_Total_NUM=(Dim_End-Dim_Begin)/Dim_Interval+1;

% read training samples
Train_DAT=zeros(NN,Train_NUM);
s=1;
for r=1:Class_NUM
   for t=1:Class_Train_NUM
     string=['E:\ORL_face\orlnumtotal\s' int2str(r) '_' int2str(t)];
     A=imread(string,'bmp');
     B=im2double(A);
     Train_DAT(:,s)=B(:);
     s=s+1;
   end
end

%read testing samples
Test_DAT=zeros(NN,Test_NUM);
s=1;
for r=1:Class_NUM
   for t=Class_Train_NUM+1:Class_Sample_NUM
     string=['E:\ORL_face\orlnumtotal\s' int2str(r) '_' int2str(t)];
     A=imread(string,'bmp');
     B=im2double(A);
     Test_DAT(:,s)=B(:);
     s=s+1;
   end
end

% to center the each training sample and testing sample
% !!! Note that: Centralization have great effection when Cos distance is
% used, but it has no impact when L2 or L1 distance is used % remove the mean value.
Mean_Image=mean(Train_DAT,2);  
Train_DAT=Train_DAT-Mean_Image*ones(1,Train_NUM);
Test_DAT=Test_DAT-Mean_Image*ones(1,Test_NUM);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Find the max Disc_NUM principal component and corresponding eigenvector
% Disc_NUM is the number of principal component,
% Projection_Matrix is the number of Projection Direction
[Projection_Matrix,disc_value]=Eigenface_f(Train_DAT, Disc_NUM);   

% LLE_DP Transform:
%Project the training simples in principal component
Train_SET=Projection_Matrix'*Train_DAT; % size of (Disc_NUM,Train_NUM); % PCA-based 

Test_SET=Projection_Matrix'*Test_DAT;   % size of (Disc_NUM,Test_NUM); % PCA-based


% for s=1:Train_NUM
%     Train_SET(:,s)=Train_SET(:,s)/norm(Train_SET(:,s));
% end
% for s=1:Test_NUM
%     Test_SET(:,s)=Test_SET(:,s)/norm(Test_SET(:,s));
% end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Train_SET=reshape(Train_SET,[Disc_NUM,Class_Train_NUM,Class_NUM]);
%Test_SET=reshape(Test_SET,[Disc_NUM,Class_Test_NUM,Class_NUM]);

A=Train_SET;
% normalization of training set
for count=1:1:200
    A(:,count)=A(:,count)/norm(A(:,count),2);
end

% Record vectors, classes the simples belong to.
classification=zeros(200,1);

% The wrong classified classes numbers.
miss=0;
% do classification work

for t=1:1:200
    y=Test_SET(:,t);
    
    % The pseudo-inverse of A
    x0=inv(A'*A)*A'*y;
    
    %solve the L1_minimization through y-Ax<=ε
    xp=l1eq_pd(x0,A,[],y,1e-3);
    %xp
    % L=norm(xp,1);% 1 Norm of xp .
    
    % 2 Norm , Record the k residual
    test=zeros(40,1);
    
    for k=1:1:40
        % Rx is the Solution of the vector, which is projected by eigenfunction.
        % The coefficient of relevant classes is non-zero. Others will be zero.
        Rx=zeros(200,1);

        % Copy the non-zero value in k class from xp to Rx. 
        % The other positions in Rx is zero because it is not copied. 
        % This equal to mapped by eigenfunction δ.
        Rx(5*(k-1)+1:5*k)=xp(5*(k-1)+1:5*k);
        
        % res=residual(y)
        res=y-A*Rx;
        
        % The 2 norm of residual in k class.
        test(k)=norm(res,2);
    end
    [value,order]=min(test);
    classification(t)=order;
   if t<(order-1)*5+1|t>order*5
       miss=miss+1;
   end
end
classification=reshape(classification,5,40);
classification
miss
Recognition_rate=(Test_NUM-miss)/Test_NUM;
Recognition_rate



%Distance_mark           :   ['Euclidean', 'L2'| 'L1' | 'Cos']
%Distance_mark='L2';
%K=3;

%for Select_DIM=Dim_Begin:Dim_Interval:Dim_End  % N is the total number of features
  
%  Miss_NUM=Classifier_NN_f(Train_SET(1:Select_DIM,:,:),Test_SET(1:Select_DIM,:,:),Distance_mark);

%  Recognition_rate=(Test_NUM-Miss_NUM)/Test_NUM;
%  fprintf('%7.3f',Recognition_rate);

%end

%fprintf('\n');

etime(clock,t0)