pca.m

使用MATLAB开发环境

clear all
close all
T=[];%所有训练图像
for i=1:40
    for j=1:5
        if (i-1)*10+j<10  %读入40*5幅训练图像
            im=imread(strcat('F:\ORL\ORL00',num2str((i-1)*10+j),'.BMP'));
        elseif (i-1)*10+j>=10 & (i-1)*10+j<100
            im=imread(strcat('F:\ORL\ORL0',num2str((i-1)*10+j),'.BMP'));
        else
            im=imread(strcat('F:\ORL\ORL',num2str((i-1)*10+j),'.BMP'));
        end
 
        temp=reshape(im,1,112*92); % 将二维图像矩阵转化为一维向量，temp是行矢量 1×N，其中N＝10304
        temp=double(temp);
        T=[T;temp];  % T是一个M*N 矩阵，T中每一行数据代表一幅图片，其中M＝200
  end
end
Tmean=mean(T); % 平均图片，1×N
for i=1:200 
    xmean(i,:)=T(i,:)-Tmean;  %xmean是一个M×N矩阵，xmean每一行保存的数据是“每幅图片数据-平均图片”
end;

sigma=xmean*xmean';   %xmean的协方差矩阵：M*M阶
[v d]=eig(sigma);     %协方差矩阵的特征向量和特征值
d1=diag(d);           %生成特征值的对角阵
[d2 index]=sort(d1);  %以升序排序
cols=size(v,2);       %特征向量矩阵的列数
for i=1:cols          %完成降序排列
    vsort(:,i)=v(:,index(cols-i+1)); % vsort 是一个M*col(注:col一般等于M)阶矩阵，保存的是按降序排列的特征向量,每一列构成一个特征向量
    dsort(i)=d1(index(cols-i+1));  % dsort保存的是按降序排列的特征值，是一维行向量
end 
%以下选择90%的能量
dsum=sum(dsort);
extractdsum=0;
p=0;
while( extractdsum/dsum < 0.9)
    p=p+1;
    extractdsum=sum(dsort(1:p));
end

%求出训练样本的特征脸，并用特征脸构成特征空间后，就可以进行人脸识别的工作了，分为训练阶段和识别阶段
% (训练阶段)计算特征脸形成的坐标系
i=1;
%选择p个较大特征值对应的特征向量
while (i<=p && dsort(i)>0)
    base(:,i)=xmean'*vsort(:,i)*dsort(i)^(-1/2);   % base是N×p阶矩阵，除以dsort(i)^(1/2)是对人脸图像的标准化
    i=i+1;
end

%将训练样本对坐标系上进行投影,即映射到特征脸形成的子空间，得到一个 M*p 阶矩阵projectedimg
projectedimg = T*base;
accu = 0;

%识别阶段
for i=1:40
    for j=6:10    %读入40x5副测试图像
        if (i-1)*10+j<10
            im=imread(strcat('F:\ORL\ORL00',num2str((i-1)*10+j),'.BMP'));
        elseif (i-1)*10+j>=10 & (i-1)*10+j<100
            im=imread(strcat('F:\ORL\ORL0',num2str((i-1)*10+j),'.BMP'));
        else
            im=imread(strcat('F:\ORL\ORL',num2str((i-1)*10+j),'.BMP'));
        end
        temp2=reshape(im,1,112*92); % 将二维图像矩阵转化为一维向量，temp2是行矢量 1×N，其中N＝10304
        temp2=double(temp2);
        projectedteimg= temp2*base; %将每一个测试样本对坐标系进行投影，是1×p阶矩阵
        for k=1:200 
            Euc_dist(k)=norm(projectedteimg-projectedimg(k,:));
        end;
        %使用三阶近邻法判别人脸图像归类 
        [Euc_dist,index2]=sort(Euc_dist);
        class1=floor(index2(1)/5)+1; %取整
        class2=floor(index2(2)/5)+1;
        class3=floor(index2(3)/5)+1;
        if class1~=class2 && class2~=class3
            class=class1;
        elseif class1==class2
            class=class1;
        elseif class2==class3
            class=class2;
        end;
        if class==i
            accu=accu+1;
           % disp(sprintf('ORL%d属于第%d组人脸图像：',(i-1)*10+j,i));
        end;
    end;
end;
disp('识别率为：');
accuracy=accu/200

%重构人脸
KL=xmean*base;
R=base*KL';
for i=1:200
    R(:,i)=R(:,i)+Tmean';
    temp3=reshape(R(:,i),112,92); 
    imwrite(mat2gray(temp3),strcat('F:\ORLRECONS\',num2str(i),'.bmp'));
end