基于PCA的人脸识别算法
--Matlab
Face recognition
Based on PCA
目录
人脸识别技术是基于人的脸部特征，对输入的人脸图象或者视频流 . 首先判断其是否存在人脸 , 如果存在人脸，则进一步的给出每个脸的位置、大小和各个主要面部器官的位置信息。并依据这些信息，进一步提取每个人脸中所蕴涵的身份特征，并将其与已知的人脸进行对比，从而识别每个人脸的身份。
广义的人脸识别实际包括构建人脸识别系统的一系列相关技术，包括人脸图像采集、人脸定位、人脸识别预处理、身份确认以及身份查找等；而狭义的人脸识别特指通过人脸进行身份确认或者身份查找的技术或系统。
主成分分析 ( Principal Component Analysis ， PCA )或者主元分析。是一种掌握事物主要矛盾的统计分析方法，它可以从多元事物中解析出主要影响因素，揭示事物的本质，简化复杂的问题。计算主成分的目的是将高维数据投影到较低维空间。给定 n 个变量的 m 个观察值，形成一个 n′m 的数据矩阵，n 通常比较大。对于一个由多个变量描述的复杂事物，人们难以认识，那么是否可以抓住事物主要方面进行重点分析呢？如果事物的主要方面刚好体现在几个主要变量上，我们只需要将这几个变量分离出来，进行详细分析。但是，在一般情况下，并不能直接找出这样的关键变量。这时我们可以用原有变量的线性组合来表示事物的主要方面， PCA 就是这样一种分析方法。
根据Kyungnam Kim在其论文 Face Recognition using Principle Component Analysis 中所陈述，PCA人脸识别的基本思想是“…express the large 1-D vector of pixelsconstructed from 2-D facial image into the compact principal components of thefeature space(eigenspace projection).” 也就是从人脸图像中找出最能代表人脸的特征空间。一个单个的人脸图片映射到这个特征空间得到这个特征空间的一组系数(这张人脸图片的特征脸特征)。如果两张人脸图片映射到这个特征空间的系数差不多，就表示这两张人脸是同一个人。从论文中我们也知道了特征空间的计算方法“Eigenspace is calculated byidentifying the eigenvectors of the covariance matrix derived from a set offacial images(vectors). ”
测试数据库为ORL_92x112，其中含40个人，每人10张图片，共400张图片。
下图为40个人的第一张图片。
6.1训练阶段
6.1.1 一张人脸图片在计算机表示为一个像素矩阵，即是一个二维数组，现在把这个二维数组变成一维数组，即把第一行后面的数全部添加到第一行。这样一张图片就能表示为一个向量d=(x1,x2......xn)。xn表示像素。
6.1.2 现在训练库里有m张人脸图片，把这些图片都表示成上述的向量形式，即d1，d2，。。。dm，把这m个向量取平均值得向量avg=(y1，y2......yn)。
得到平均图像:
6.1.3 用d1，d2...........dm分别减去avg后组成一个矩阵A,即矩阵A的第一行为d1-avg，后面类似。A的大小为m×n。因为找特征空间不能基于一张图片，而要在所有的人脸图片提取出共同特征，所以要取各个人脸向量到平均人脸向量的向量差。依据这个每个人脸图片跟平均脸向量的向量差组成矩阵A，然后依据矩阵A来求解最特征空间。
6.1.4.矩阵A乘以A的逆矩阵A‘得A的协方差矩阵B，B的大小m×m，求B的特征向量。取最大的K个特征向量组成新的矩阵T，T的大小m×k。
6.1.5 使用A’乘以T得到特征脸C，C的大小n×k。
6.1.6 用图片向量d乘以C得到图片向量d在特征脸的投影向量pn，有多少张图片就有多少个pn。pn的大小1×k
6.2探究各特征值所占有的能量数
6.2.1 第n个特征值占所有特征值之和的百分比 (特征值从小到大排列)。
6.2.2前n个特征值占所有特征值之和的百分比 (特征值从小到大排列)。
在实验中，要求保留90%的能量，所以当取前57个特征值得时候就可以满足要求了。
6.3测试阶段
6.3.1.一张新的图片也表示为d的向量，记为D，D的大小1×n
6.3.2. D乘以上面训练得到的特征脸C得到这个图片向量D在C下的投影向量P，p的大小1×k。
6.3.3.计算p与上面所有的pn的向量距离，与p最小的那个向量所对应的人脸图片跟这张新人脸图片最像。
6.3.4.判别图像的相似性有两种方法：一种是计算N维空间中图像间的距离，另一种是测量图像间的相似性。当测量距离时，距离应尽可能的小，一般选择距离测试图像最近的训练图像作为它所属的类别，识别的方法和最初的图像匹配方法类似：将测试集中的每一幅降维图像与降维的训练集进行匹配，然后将其分类到距离最小的训练集头像中，如果两个头像表示一个人，表示识别成功。而测量相似性时，图像应尽可能的相似，也就是说具有最大相似性的训练图像类被认为是测试图像所属的类别。本程序采用三阶近邻法。核心代码如下
运行程序后，将读取40个人的后5张图片，并进行比对，代码附录给出。
可得到识别率为accuracy =0.8500；
[1]Kyungnam Kim. Face Recognition using Principle Component Analysis [J] . IEEESignal Processing Society，2002，9(2)：40-42.
[2] 罗韵. 用Matlab进行PCA人脸识别[N/OL] . 网易LOFTER，[-04-12] ./luoyun_dreamer/blog/static/21552907031281538309/?suggestedreading
[3]haitao111313 . 基于PCA的人脸检测 [N/OL] . CSDN博客 , [-08-16] ./haitao111313/article/details/7875392
[4] Turk,M.A. ; Media Lab., MIT, Cambridge, MA, USA ; Pentland, A.P. Face Recognitionusing Eigenfaces [J] . IEEE Signal Processing Society，1991，7
附录1.基于PCA的人脸识别的Matlab实现代码
function FaceRecognition
clear % calc xmean,sigma and its eigen decomposition
close all
allsamples=[];%所有训练图像
for i=1:40
for j=1:5
if(i<10)
a=imread(strcat('E:\ORL_92x112\00',num2str(i),'0',num2str(j),'.bmp'));
else
a=imread(strcat('E:\ORL_92x112\0',num2str(i),'0',num2str(j),'.bmp'));
end
b=a(1:112*92); % b是行矢量 1×N，其中N＝10304，提取顺序是先列后行，即从上到下，从左到右
b=double(b);
allsamples=[allsamples; b]; % allsamples 是一个M * N 矩阵，allsamples 中每一行数据代表一张图片，其中M＝200
end
end
samplemean=mean(allsamples); % 平均图片，1 × N
figure%平均图
imshow(mat2gray(reshape(samplemean,112,92)));
for i=1:200
xmean(i,:)=allsamples(i,:)-samplemean; % xmean是一个M × N矩阵，xmean每一行保存的数据是“每个图片数据-平均图片”
end;
% figure%平均图
% imshow(mat2gray(reshape(xmean(1,:),112,92)));
sigma=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);
dsum_extract = 0;
p = 0;
while( dsum_extract/dsum < 0.9)
p = p + 1;
dsum_extract = sum(dsort(1:p));
end
a=1:1:200;
for i=1:1:200
y(i)=sum(dsort(a(1:i)) );
end
figure
y1=ones(1,200);
plot(a,y/dsum,a,y1*0.9,'linewidth',2);
grid
title('前n个特征特占总的能量百分比');
xlabel('前n个特征值');
ylabel('占百分比');
figure
plot(a,dsort/dsum,'linewidth',2);
grid
title('第n个特征特占总的能量百分比');
xlabel('第n个特征值');
ylabel('占百分比');
i=1; % (训练阶段)计算特征脸形成的坐标系
while (i<=p && dsort(i)>0)
base(:,i) = dsort(i)^(-1/2) * xmean' * vsort(:,i); % base是N×p阶矩阵，除以dsort(i)^(1/2)是对人脸图像的标准化，特征脸
i = i + 1;
end % add by wolfsky 就是下面两行代码，将训练样本对坐标系上进行投影,得到一个 M*p 阶矩阵allcoor
%
% for i=1:20
% figure%平均图
% b=reshape(base(:,i)',112,92);%
% imshow(mat2gray(b));
% end
allcoor = allsamples * base; accu = 0; % 测试过程
for i=1:40
for j=6:10 %读入40 x 5 副测试图像
if(i<10)
if(j<10)
a=imread(strcat('E:\ORL_92x112\00',num2str(i),'0',num2str(j),'.bmp'));
else
a=imread(strcat('E:\ORL_92x112\00',num2str(i),num2str(j),'.bmp'));
end
else
if(j<10)
a=imread(strcat('E:\ORL_92x112\0',num2str(i),'0',num2str(j),'.bmp'));
else
a=imread(strcat('E:\ORL_92x112\0',num2str(i),num2str(j),'.bmp'));
end
end
b=a(1:10304);
b=double(b);
tcoor= b * base; %计算坐标，是1×p阶矩阵
for k=1:200
mdist(k)=norm(tcoor-allcoor(k,:));
end; %三阶近邻
[dist,index2]=sort(mdist);
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;
end;
end;
end;
accuracy=accu/200 %输出识别率
