my_2dpca.m

现有的代数特征的抽取方法绝大多数采用一维的方法

function recognize = my_2dpca(path,persons,numface,numtrain,numpc)
% RECOGNIZE = MY_2DPCA(PATH,PERSONS,NUMFACE,NUMTRAIN,NUMPC)
% Recognize 正确识别率; path 存放人脸图像的路径; 
% persons 类别数; numface 每个类别的样本数;numtrain 训练样本的个数; 
% numpc 主元个数

% 例：用orl人脸库，训练样本为8，选取主元数目为100
% path='D:\orl\';
% persons=40;
% numface=10;
% numtrain=6;
% numpc=6;
% recognize=my_2dpca(path,persons,numface,numtrain,numpc)
% 结果显示：    正确识别率：96.5%

% ********************************************************************************************************
% Part I 环境设置
% ********************************************************************************************************

% 《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 >>> 计算参数   ====================================================================================

    numtest=numface-numtrain;                                    % 每个人参加测试的样本数
    all_train=persons*numtrain;                                  % 所有参加训练的样本数 
    all_test=persons*numtest;                                    % 所有参加测试的样本数
    
% ********************************************************************************************************
% Part II 求正交坐标系(由正交的单位向量{X1,X2,...,Xnumpc}构成的空间坐标系)
% ********************************************************************************************************

% <<< 1 >>> 装载训练样本 A(:,:,i) ========================================================================

    for i=1:persons                                          
        for j=1:numtrain                                        
            img_path=[Folder(i,:),BMPfile(j,:)];           
	        img=double(imread(img_path));                   % 读入图像并转换为双精度以便计算
            if i==1 & j==1                                  % 用第一个样本来确定人脸样本的大小   
                Dim=size(img);
            end
            A(:,:,(i-1)*numtrain+j)=img;                    % A 的第三维用作训练样本的序号
        end
    end

% <<< 2 >>> 计算平均脸 ef 及训练样本规范化 ===============================================================
    
    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 >>> 计算协方差矩阵 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 为特征值构成的对角阵 ( col x col )

% <<< 6 >>> 找出最大(主元)的 numpc 个特征值对应的特征向量(相互正交的单位向量)作基构成坐标系 =================

    [ld,ind]=sort(diag(D));                         % 对特征值进行索引排序
    landa=ld(end:-1:1);                             % 从大到小排序
    Ind=ind(end:-1:1);                              % 把对应的序号排序
    for i=1:numpc
        X(:,i)=V(:,Ind(i));                         % 由大到小排前 numpc 个特征值的特征向量
    end


% ********************************************************************************************************
% Part III 求 A , B 的坐标(向量组)
% ********************************************************************************************************

% <<< 1 >>> 每个训练样本A(:,:,i)投影到正交的单位向量{X1,X2,...,Xnumpc}构成的空间坐标上,
%           得到它的坐标---投影特征向量组{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  在坐标系内的欧氏距离,并进行分类
% ********************************************************************************************************

err=0;                                      % 错判计数器
Aclass=zeros(all_train,1);                        % 训练样本的分类
Bclass=zeros(all_test,1);                        % 测试样本的本应属于的类别
Bresult=zeros(all_test,1);                       % 测试样本被判别属于的类别

for ac=1:all_train
    Aclass(ac)=ceil(ac/numtrain);                % 训练样本的分类
end

for bc=1:all_test
    Bclass(bc)=ceil(bc/numtest);                % 测试样本的本应类别
end

for i=1:all_test
    Dij=zeros(all_train,1);                       % 第 i 个测试样本到所有训练样本的距离
    for j=1:all_train
        YD=YT(:,:,i)-Y(:,:,j);              % 第 i 个测试样本与第 j 个训练样本的各个投影特征向量的差
        for k=1:numpc
            Dij(j)=Dij(j)+norm(YD(:,k));    % 差的二范数累加即为两者的欧氏距离
        end
    end
    [Min,IND]=sort(Dij);                    % 从小到大索引,取最小的那个类别即为结果
    Bresult(i)=Aclass(IND(1));
    if Bresult(i)~=Bclass(i)                % 判别错误率
        err=err+1;
    end
end
    
recognize=['正确识别率:',num2str((1-err/all_test)*100),'%'];    % 显示识别结果