基于PCA的人脸特征提取和识别.doc

工程设计报告设计题目：基于PCA的人脸特征抽取与识别学 院： 电子工程学院 专 业： 智能科学与技术 ： 钟佩 学号： 02085156 时 间： 2011年11月 指导教师： 缑水平 目录摘要31PCA进行特征抽取和识别的方法与理论基础31.1 K-L 变换31.2 利用 PCA 进行人脸识别41.3 PCA 的理论基础41.3.1 投影41.3.2 PCA 的作用与其统计特性51.3.3 特征脸61.3.4 图片重建71.3.5 奇异值分解（SVD）71.3.6 利用小矩阵计算大矩阵特征向量81.3.7 图片归一化82结果91.识别率92.特征脸93.人脸重构103参考文献104 附录matlab 源码114.1 人脸识别114.2 特征人脸124.3 人脸重建14摘要对于一幅图像可以看作一个由像素值组成的矩阵，也可以扩展开，看成一个矢量，如一幅N*N 象素的图像可以视为长度为N2 的矢量，这样就认为这幅图像是位于N2 维空间中的一个点，这种图像的矢量表示就是原始的图像空间，但是这个空间仅是可以表示或者检测图像的许多个空间中的一个。不管子空间的具体形式如何，这种方法用于图像识别的基本思想都是一样的，首先选择一个合适的子空间，图像将被投影到这个子空间上，然后利用对图像的这种投影间的某种度量来确定图像间的相似度，最常见的就是各种距离度量。在本文中，我们将讨论PVA算法来对人脸进行特征抽取和识别。1PCA进行特征抽取和识别的方法与理论基础1.1 K-L 变换PCA 方法是由Turk 和Pentlad 提出来的，它的基础就是Karhunen-Loeve 变换(简称KL变换)，是一种常用的正交变换。下面我们首先对K-L 变换作一个简单介绍:假设X 为n 维的随机变量，X 可以用n 个基向量的加权和来表示:式中: i是加权系数，i是基向量，此式还可以用矩阵的形式表示:取基向量为正交向量，即则系数向量为：= T X综上所述，K-L 展开式的系数可用下列步骤求出:步骤一求随即向量X 的自相关矩阵，由于没有类别信息的样本集的均值向量，常常没有意义，所以也可以把数据的协方差矩阵作为K-L 坐标系的产生矩阵，这里是总体均值向量。步骤二求出自相关矩阵或协方差矩阵R 的本征值和本征向量，步骤三展开式系数即为= T XK_L 变换的实质是建立了一个新的坐标系，将一个物体主轴沿特征矢量对齐的旋转变换，这个变换解除了原有数据向量的各个分量之间相关性，从而有可能去掉那些带有较少信息的坐标系以达到降低特征空间维数的目的。1.2 利用 PCA 进行人脸识别完整的PCA 人脸识别的应用包括几个步骤：人脸图像预处理；读入人脸库，训练形成特征子空间；把训练图像和测试图像投影到上一步骤中得到的子空间上；选择一定的距离函数进行识别。下面详细描述整个过程（源码见faceRec.m）。1. 读入人脸库归一化人脸库后，将库中的每人选择一定数量的图像构成训练集，其余构成测试集。设归一化后的图像是n*m,按列相连就构成N=n*m 维矢量，可视为N 维空间中的一个点，可以通过K-L 变换用一个低维子空间描述这个图像。2. 计算K- L 变换的生成矩阵所有训练样本的协方差矩阵为（以下三个等价）：是平均人脸, M 训练人脸数，协方差矩阵CA是一个N*N的矩阵, N 是xi的维数。为了方便计算特征值和特征向量，一般选用第2个公式。根据K - L 变换原理，我们所求的新坐标系即由矩阵A·AT的非零特征值所对应的特征向量组成。直接求N*N大小矩阵CA的特征值和正交归一特征向量是很困难的, 根据奇异值分解原理（见段落1.2.5和1.2.6），可以通过求解AT·A的特征值和特征向量来获得AT·A的特征值和特征向量。在计算得到CA的所有非零特征值（从大到小排序，1 r < M ）与其对应的单位正交特征向量后，可以得到特征空间，从而可以计算一图片X在特征空间上的投影系数（也可以理解为X在空间U中的坐标）：3. 识别利用公式（2），首先把所有训练图片进行投影，然后对于测试图片也进行同样的投影，采用判别函数对投影系数进行识别。1.3 PCA 的理论基础1.3.1 投影设d 维样本x1，x 2，,xn，以与一个d 维基w ，那么标量：yi=wT xi是相当于xi在基上的坐标值。如果w = 1 , yi就是把i x 向方向为w 的直线进行投影的结果，可以从图1 看到。推广之，如果有一组基（m 个）组成的空间，那么可以得到xi在空间W上的坐标为：。图1 投影图进一步，表达式w = m+ ae表示w 是一条通过点m，方向为e的直线。1.3.2 PCA 的作用与其统计特性采用PCA 对原始数据的处理，通常有三个方面的作用降维、相关性去除、概率估计。下面分别进行介绍：去除原始数据相关性从统计学上讲，EX E(X )Y E(Y) 称为随机变量X 与Y 协方差，记为Cov(X,Y)。令，称为随机变量X 与Y 的相关系数。XY= 1则X 与Y 是相关的，XY= 0，则X 与Y 是不相关的。命题1 对于矩阵A 来说，如果AAT是一个对角阵，那么A中的向量是非相关的。由PCA 处理的人脸库数据的非相关性可以从两点进行说明。（1）基底的非相关性特征空间基是非相关的，即UUT= I 。（2）投影系数的非相关性由SVD 可知, 其中, mx是平均人脸。根据公式（2）可以把A 映射到特征空间上，得到：B =UT * A，其中B 是非相关的，可由下面得到证明：Y的协方差矩阵为：由命题1 可知，B 是非相关的。统计参数（均值与方差）均值即mx-平均人脸。命题2 随机变量方差越大，包含的信息越多，当一个变量方差为0 时，该变量为常数，不含任何信息。用PCA 计算主分量，就是寻找一组向量，使得原始数据A = 1,2 ,.,M 在这组向量上的投影值的方差尽可能大。最大方差对应的向量就是第一主成份，以后递推就是第二主成份，第三主成份。用PCA 计算主分量就是求原始数据A = 1,2 ,.,M （其中i= ximx）协方差矩阵的特征向量，由公式（3）可知，是A在ui上的投影值，其中P 的方差就ui对应的特征值i，可以理解为：命题3 所有原始数据在主分量ui上的投影值方差为i。降维如果在原始空间表示一幅n*m 大小的图片X，那么需要一个Nn*m 维矢量，但是当用公式(2)把它映射到特征空间后，只需要一个r*1 维的向量就可。另外，由命题2 可知，可以根据方差的大小来判断特征向量的重要性。由ORL 图片库的200 个人脸计算得到的特征值呈图2 分布，可知特征向量重要性呈指数下降，据此可以只选用前面几个重要的特征向量来构建特征空间。通过计算，前71 个特征值占了90.17%，因此r 可以取71 而非200，从而达到进一步降维的作用。图2 特征值的分布1.3.3 特征脸中的每一个单位向量都构成一个特征脸，如图3 所示。由这些特征脸所成的空间称为特征脸子空间，需要注意对于正交基的选择的不同考虑，对应较大特征值的特征向量(正交基)也称主分量，用于表示人脸的大体形状，而对应于较小特征值的特征向量则用于描述人脸的具体细节，或者从频域来看，主分量表示了人脸的低频部分，而此分量则描述了人脸的高频部分（源码见EigenFace.m）。 1 2 10 50 70 平均脸图3 特征脸，分别是第1，2，10，50，70 分量，最后一是平均脸。1.3.4 图片重建要进行图片X的重建，首先对X投影到特征空间上，得到系数,然后选用一部分系数与特征向量进行原始图片的重建：，其中1: t表示取前t 个元素。（见reconstruct.m）在图4 中，其中前两图片来自训练样本，第3 来自测试样本，可以看到对于训练样本，PCA 系数可以对图片实现很好重建，而对于训练样本以外的图片重建效果很差。原始图像 t=15 t=50 t=100 t=150 t=199图4 人脸图像重建。第1列图片是输入原始图，其它列图片是重建结果，数字表示t 的数目。1.3.5 奇异值分解（SVD）设A 是秩为r 的m*n（m>>n）维矩阵，则存在两个正交矩阵和一个对角阵：其中，且UUT= I ，VVT= I ，i呈降序排列。其中i2为AAT和ATA 的非零特征值，ui和vi分别是AAT和ATA对应于i2的特征向量。可得一个推论：可以计算ATA的特征值i2与相应的正交归一特征向量i v 后，可由推论知AAT 的正交归一特征向量 注意，协方差矩阵的特征值为：i2/M 。1.3.6 利用小矩阵计算大矩阵特征向量高阶矩阵的特征向量可以转化为求低阶矩阵的特征向量：设：A 是秩为r 的m*n（m>>n）维矩阵，是一个矩阵，现在要求Cx的特征值与特征向量，可通过先求小矩阵的特征向量和特征值，两者之间有以下关系：显然,Cx= AAT的特征向量是Avi（注意没有单位化），亦为其特征值。结论：1.2.5 与1.2.6 的方法计算协方差矩阵的特征向量，特征值的结果是一致的，只是要注意1.2.5 中的特征值要除以M，1.2.6 中的特征向量要单位化。1.3.7 图片归一化图片标准化通常是一个整体概念，要求把图片归一到均值为0，方差为1 下情况下。这个概念类似于一般正态分布向标准正态分布的转化：命题4 若X N(,2)，则所以要对一组图片中的一Xi进行归一化（标准化），只需要减去均值，除以方差就可以了。均值，方差为2结果1.识别率2.特征脸3.人脸重构3参考文献1 数字图像处理与机器视觉Visual C+与Matlab实现 铮 王艳平 薛桂香著4附录matlab 源码4.1 人脸识别% FaceRec.m% PCA 人脸识别修订版，识别率88% calc xmean,sigma and its eigen decompositionallsamples=;%所有训练图像for i=1:40for j=1:5a=imread(strcat('D:ORLs',num2str(i),'',num2str(j),'.pgm');% imshow(a);b=a(1:112*92); %b是行矢量 1×N，其中N10304，提取顺序是先列后行，即从上到下，从左到右b=double(b);allsamples=allsamples; b; %allsamples 是一个M * N 矩阵，allsamples 中每一行数据代表一图片，其中M200endendsamplemean=mean(allsamples); % 平均图片，1 ×Nfor i=1:200 xmean(i,:)=allsamples(i,:)-samplemean; % xmean 是一个M ×N 矩阵，xmean每一行保存的数据是“每个图片数据-平均图片”end;% 获取特征值与特征向量sigma=xmean*xmean' % M * M 阶矩阵v d=eig(sigma);d1=diag(d);% 按特征值大小以降序排列dsort = flipud(d1);vsort = fliplr(v);%以下选择90%的能量dsum = sum(dsort);dsum_extract = 0;p = 0;while( dsum_extract/dsum < 0.9)p = p + 1;dsum_extract = sum(dsort(1:p);endi=1;%(训练阶段)计算特征脸形成的坐标系base = xmean' * vsort(:,1:p) * diag(dsort(1:p).(-1/2);%base 是N×p 阶矩阵，除以dsort(i)(1/2)是对人脸图像的标准化（使其方差为1）% xmean' * vsort(:,i)是小矩阵的特征向量向大矩阵特征向量转换的过程%while (i<=p && dsort(i)>0)%base(:,i) = dsort(i)(-1/2) * xmean' * vsort(:,i); % base 是N×p 阶矩阵，除以dsort(i)(1/2)是对人脸图像的标准化（使其方差为1）% i = i + 1; %xmean' * vsort(:,i)是小矩阵的特征向量向大矩阵特征向量转换的过程%end% 以下两行add by gongxun 将训练样本对坐标系上进行投影,得到一个M*p 阶矩阵allcoorallcoor = allsamples * base; % allcoor 里面是每训练人脸图片在M*p 子空间中的一个点，即在子空间中的组合系数accu = 0; %下面的人脸识别过程中就是利用这些组合系数来进行识别% 测试过程for i=1:40for j=6:10 %读入40 x 5 副测试图像a=imread(strcat('D:ORLs',num2str(i),'',num2str(j),'.pgm');b=a(1:10304);b=double(b);tcoor= b * base; %计算坐标，是1×p 阶矩阵for k=1:200mdist(k)=norm(tcoor-allcoor(k,:);end;%三阶近邻dist,index2=sort(mdist);class1=floor( (index2(1)-1)/5 )+1;class2=floor(index2(2)-1)/5)+1;class3=floor(index2(3)-1)/5)+1;if class1=class2 && class2=class3class=class1;elseif class1=class2class=class1;elseif class2=class3class=class2;end;if class=iaccu=accu+1;end;end;end;accuracy=accu/200 %输出识别率4.2 特征人脸% eigface.mfunction = eigface()% calc xmean,sigma and its eigen decompositionclose all;clc;clear all;allsamples=;%所有训练图像for i=1:40for j=1:5a=imread(strcat('D:ORLs',num2str(i),'',num2str(j),'.pgm');% imshow(a);b=a(1:112*92); % b 是行矢量1×N，其中N10304，提取顺序是先列后行，即从上到下，从左到右b=double(b);allsamples=allsamples; b; %allsamples 是一个M * N 矩阵，allsamples 中每一行数据代表一图片，其中M200endendsamplemean=mean(allsamples); % 平均图片，1 ×Nfor i=1:200 xmean(i,:)=allsamples(i,:)-samplemean; % xmean 是一个M ×N 矩阵，xmean每一行保存的数据是“每个图片数据-平均图片”end;% 获取特征值与特征向量sigma=xmean*xmean' % M * M 阶矩阵v d=eig(sigma);d1=diag(d);%按特征值大小以降序排列dsort = flipud(d1);vsort = fliplr(v);%以下选择90%的能量dsum = sum(dsort);dsum_extract = 0;p = 0;while( dsum_extract/dsum < 0.9)p = p + 1;dsum_extract = sum(dsort(1:p);endp = 199;% (训练阶段)计算特征脸形成的坐标系%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; % xmean' * vsort(:,i)是小矩阵的特征向量向大矩阵特征向量转换的过程%endbase = xmean' * vsort(:,1:p) * diag(dsort(1:p).(-1/2);% 生成特征脸,并显示出来for (k=1:p),temp = reshape(base(:,k), 112,92);if(k<21)figure(1);subplot(4,5,k),imshow(mat2gray(temp);elseif(k<41) figure(2); subplot(4,5,(k-20),imshow(mat2gray(temp);elseif(k<61) figure(3); subplot(4,5,(k-40),imshow(mat2gray(temp);elseif(k<81) figure(4); subplot(4,5,(k-60),imshow(mat2gray(temp);elseif(k<101) figure(5); subplot(4,5,(k-80),imshow(mat2gray(temp); elseif(k<121) figure(6); subplot(4,5,(k-100),imshow(mat2gray(temp); elseif(k<141) figure(7); subplot(4,5,(k-120),imshow(mat2gray(temp);elseif(k<161) figure(8); subplot(4,5,(k-140),imshow(mat2gray(temp); elseif(k<181) figure(9); subplot(4,5,(k-160),imshow(mat2gray(temp);else(k<200) figure(10); subplot(4,5,(k-180),imshow(mat2gray(temp);endendavg = reshape(samplemean, 112,92);subplot(4,5,20),imshow(mat2gray(avg),title('平均脸');%将模型保存save('D:ORLmodel.mat', 'base', 'samplemean');4.3 人脸重建% Reconstruct.mfunction = reconstruct()close all;clear all;clcload D:ORLmodel.mat;%计算新图片在特征子空间中的系数img = 'D:ORLs23.pgm'a=imread(img);b=a(1:112*92); % b 是行矢量1×N，其中N10304，提取顺序是先列后行，即从上到下，从左到右b=double(b);b=b-samplemean;c = b * base; % c 是图片a 在子空间中的系数, 是1*p 行矢量% 根据特征系数与特征脸重建图figure(1)subplot(1,6,1),imshow(img);% 前15 个t = 15;temp = base(:,1:t) * c(1:t)'temp = temp + samplemean'subplot(1,6,2),imshow(mat2gray(reshape(temp, 112,92);% 前50 个t = 50;temp = base(:,1:t) * c(1:t)'temp = temp + samplemean'subplot(1,6,3),imshow(mat2gray(reshape(temp, 112,92);% 前100 个t = 100;temp = base(:,1:t) * c(1:t)'temp = temp + samplemean'subplot(1,6,4),imshow(mat2gray(reshape(temp, 112,92);% 前150 个t = 150;temp = base(:,1:t) * c(1:t)'temp = temp + samplemean'subplot(1,6,5),imshow(mat2gray(reshape(temp, 112,92);% 前199 个t = 199;temp = base(:,1:t) * c(1:t)'temp = temp + samplemean'subplot(1,6,6),imshow(mat2gray(reshape(temp, 112,92);14 / 14