function recognize = my_2dpca(path,persons,numface,numtrain,numpc)
% RECOGNIZE = MY_2DPCA(PATH,PERSONS,NUMFACE,NUMTRAIN,NUMPC)
% Recognize 正確識別率; path 存放人臉圖像的路徑; 
% persons 類別數(shù); numface 每個類別的樣本數(shù);numtrain 訓(xùn)練樣本的個數(shù); 
% numpc 主元個數(shù)

% 例：用orl人臉庫，訓(xùn)練樣本為8，選取主元數(shù)目為100
% path='D:\orl\';
% persons=40;
% numface=10;
% numtrain=6;
% numpc=6;
% recognize=my_2dpca(path,persons,numface,numtrain,numpc)
% 結(jié)果顯示：    正確識別率：96.5%

% ********************************************************************************************************
% Part I 環(huán)境設(shè)置
% ********************************************************************************************************

% 《1》生成路徑 ==========================================================================================

for i=1:persons
    str=num2str(i);
    if i<10
        str=['0',str];
    end
    Folder(i,:)=[path,'s',str,'\'];
end
 
BMPfile =['01.bmp';'02.bmp';'03.bmp';'04.bmp';'05.bmp';
     '06.bmp';'07.bmp';'08.bmp';'09.bmp';'10.bmp';'11.bmp'];

% <<< 2 >>> 計算參數(shù)   ====================================================================================

    numtest=numface-numtrain;                                    % 每個人參加測試的樣本數(shù)
    all_train=persons*numtrain;                                  % 所有參加訓(xùn)練的樣本數(shù) 
    all_test=persons*numtest;                                    % 所有參加測試的樣本數(shù)
    
% ********************************************************************************************************
% Part II 求正交坐標系(由正交的單位向量{X1,X2,...,Xnumpc}構(gòu)成的空間坐標系)
% ********************************************************************************************************

% <<< 1 >>> 裝載訓(xùn)練樣本 A(:,:,i) ========================================================================

    for i=1:persons                                          
        for j=1:numtrain                                        
            img_path=[Folder(i,:),BMPfile(j,:)];           
	        img=double(imread(img_path));                   % 讀入圖像并轉(zhuǎn)換為雙精度以便計算
            if i==1 & j==1                                  % 用第一個樣本來確定人臉樣本的大小   
                Dim=size(img);
            end
            A(:,:,(i-1)*numtrain+j)=img;                    % A 的第三維用作訓(xùn)練樣本的序號
        end
    end

% <<< 2 >>> 計算平均臉 ef 及訓(xùn)練樣本規(guī)范化 ===============================================================
    
    ef=zeros(Dim(1),Dim(2));
    for i=1:all_train
        ef=ef+A(:,:,i);                         
    end
    ef=ef/all_train;                                 
    for i=1:all_train
        A(:,:,i)=A(:,:,i)-ef;                  
    end
       
   
% <<< 3 >>> 計算協(xié)方差矩陣 G ============================================================================

    G=zeros(Dim(2),Dim(2));    
    for i=1:all_train
        G=G+A(:,:,i)'*A(:,:,i);                                   
    end
    G=G/all_train;
    
% <<< 5 >>> 計算特征值 D ,特征向量 V ====================================================================

    [V,D]=eig(G);                                   % V 為特征向量矩陣 ( col x col ) 
                                                    % D 為特征值構(gòu)成的對角陣 ( col x col )

% <<< 6 >>> 找出最大(主元)的 numpc 個特征值對應(yīng)的特征向量(相互正交的單位向量)作基構(gòu)成坐標系 =================

    [ld,ind]=sort(diag(D));                         % 對特征值進行索引排序
    landa=ld(end:-1:1);                             % 從大到小排序
    Ind=ind(end:-1:1);                              % 把對應(yīng)的序號排序
    for i=1:numpc
        X(:,i)=V(:,Ind(i));                         % 由大到小排前 numpc 個特征值的特征向量
    end


% ********************************************************************************************************
% Part III 求 A , B 的坐標(向量組)
% ********************************************************************************************************

% <<< 1 >>> 每個訓(xùn)練樣本A(:,:,i)投影到正交的單位向量{X1,X2,...,Xnumpc}構(gòu)成的空間坐標上,
%           得到它的坐標---投影特征向量組{Y1(i),Y2(i),...Ynumpc(i)}

    for i=1:all_train    
        Y(:,:,i)=A(:,:,i)*X;                        % Y 為 ( row x numpc ) 矩陣
    end

% <<< 2 >>> 裝載測試樣本 B(:,:,i) , 減去平均臉 ef , 計算在 X 空間的坐標 ================================

    for i=1:persons                                          
        for j=(numtrain+1):numface                      % 從每個人第(numtrain+1)個樣本起剩余的樣本用來測試                              
            img_path=[Folder(i,:),BMPfile(j,:)];           
	        img=double(imread(img_path));           
            B(:,:,numtest*(i-1)+(j-numtrain))=img;          
        end
    end
    for i=1:all_test
        B(:,:,i)=B(:,:,i)-ef;                       % 減去平均臉
    end
    for i=1:all_test    
        YT(:,:,i)=B(:,:,i)*X;                       % YT 為 B 在坐標系 {X1,X2,...,Xnumpc} 上的坐標
    end

% ********************************************************************************************************
% Part IV 計算的 A , B  在坐標系內(nèi)的歐氏距離,并進行分類
% ********************************************************************************************************

err=0;                                      % 錯判計數(shù)器
Aclass=zeros(all_train,1);                        % 訓(xùn)練樣本的分類
Bclass=zeros(all_test,1);                        % 測試樣本的本應(yīng)屬于的類別
Bresult=zeros(all_test,1);                       % 測試樣本被判別屬于的類別

for ac=1:all_train
    Aclass(ac)=ceil(ac/numtrain);                % 訓(xùn)練樣本的分類
end

for bc=1:all_test
    Bclass(bc)=ceil(bc/numtest);                % 測試樣本的本應(yīng)類別
end

for i=1:all_test
    Dij=zeros(all_train,1);                       % 第 i 個測試樣本到所有訓(xùn)練樣本的距離
    for j=1:all_train
        YD=YT(:,:,i)-Y(:,:,j);              % 第 i 個測試樣本與第 j 個訓(xùn)練樣本的各個投影特征向量的差
        for k=1:numpc
            Dij(j)=Dij(j)+norm(YD(:,k));    % 差的二范數(shù)累加即為兩者的歐氏距離
        end
    end
    [Min,IND]=sort(Dij);                    % 從小到大索引,取最小的那個類別即為結(jié)果
    Bresult(i)=Aclass(IND(1));
    if Bresult(i)~=Bclass(i)                % 判別錯誤率
        err=err+1;
    end
end
    
recognize=['正確識別率:',num2str((1-err/all_test)*100),'%'];    % 顯示識別結(jié)果