龚 勋(2007-4-29)
1.1 问题描述[1]
对于一幅图像可以看作一个由像素值组成的矩阵,也可以扩展开,看成一个矢量,如一幅N*N象素的图像可以视为长度为N2的矢量,这样就认为这幅图像是位于N2维空间中的一个点,这种图像的矢量表示就是原始的图像空间,但是这个空间仅是可以表示或者检测图像的许多个空间中的一个。不管子空间的具体形式如何,这种方法用于图像识别的基本思想都是一样的,首先选择一个合适的子空间,图像将被投影到这个子空间上,然后利用对图像的这种投影间的某种度量来确定图像间的相似度,最常见的就是各种距离度量。
1.1.1 K-L变换[1]
PCA方法是由Turk和Pentlad提出来的,它的基础就是Karhunen-Loeve变换(简称KL变换),是一种常用的正交变换。下面我们首先对K-L变换作一个简单介绍:
假设X为n维的随机变量,X可以用n个基向量的加权和来表示:
X=∑αiφi
i=1n
式中: αi是加权系数,φi是基向量,此式还可以用矩阵的形式表示:
X=(φ1,φ2,,φn)(α1,α2,,αn)T=Φα
取基向量为正交向量,即
⎧1 i=j
ΦTΦj=⎨⇒ΦTΦj=I
⎩0 i≠j
则系数向量为:
α=ΦTX
综上所述,K-L展开式的系数可用下列步骤求出:
T
步骤一 求随即向量X的自相关矩阵R=E⎡⎣XX⎤⎦,由于没有类别信息的样本集的µ均值向T
量,常常没有意义,所以也可以把数据的协方差矩阵∑=E⎡⎣(x−µ)(x−µ)⎤⎦作为K_L坐标系的产生矩阵,这里µ是总体均值向量。
步骤二 求出自相关矩阵或协方差矩阵R的本征值λi和本征向量φi,Φ=(φ1,φi,,φn)
步骤三 展开式系数即为α=ΦTX
K_L变换的实质是建立了一个新的坐标系,将一个物体主轴沿特征矢量对齐的旋转 变换,这个变换解除了原有数据向量的各个分量之间相关性,从而有可能去掉那些带 有较少信息的坐标系以达到降低特征空间维数的目的。
1.1.2 利用PCA进行人脸识别
完整的PCA人脸识别的应用包括几个步骤:人脸图像预处理;读入人脸库,训练形成特征子空间;把训练图像和测试图像投影到上一步骤中得到的子空间上;选择一定 的距离函数进行识别。下面详细描述整个过程(源码见’faceRec.m’)。 1. 读入人脸库
归一化人脸库后,将库中的每人选择一定数量的图像构成训练集,其余构成测试集。设归一化后的图像是n*m,按列相连就构成N=n*m维矢量,可视为N维空间中的一个点,可以通过K-L变换用一个低维子空间描述这个图像。 2. 计算K- L变换的生成矩阵
所有训练样本的协方差矩阵为(以下三个等价):
M
⎧
=1. C(xkixkT)/M−mximxT∑A⎪
k=1
⎪⎪T
(1) ⎨2. CA=(AiA)/M
⎪M
⎪3.CA=⎡∑(xi−mx)(xi−mx)T⎤M⎢⎥⎪i=1⎣⎦⎩
A={φ1,φ2,...,φM}, φi=xi−mx,mx是平均人脸, M 训练人脸数,协方差矩阵CA 是
一个N*N的矩阵, N 是xi的维数。
为了方便计算特征值和特征向量,一般选用第2个公式。根据K - L 变换原理,我们所求的新坐标系即由矩阵AiAT的非零特征值所对应的特征向量组成。直接求N*N大小矩阵CA的特征值和正交归一特征向量是很困难的, 根据奇异值分解原理(见段落1.2.5和1.2.6),可以通过求解ATiA的特征值和特征向量来获得ATiA的特征值和特征向量,。
在计算得到CA的所有非零特征值[λ0,λ1,,λr−1](从大到小排序,1≤r 3. 识别 利用公式(2),首先把所有训练图片进行投影,然后对于测试图片也进行同样的投影,采用判别函数对投影系数进行识别。 1.2 PCA的理论基础 1.2.1 投影[2] 设d维样本x1,x2,,xn,以及一个d维基w,那么标量: yi=wTxi 是相当于xi在基上的坐标值。如果w=1,yi就是把xi向方向为w的直线进行投影的结果,可以从图 1看到。推广之,如果有一组基(m个)组成的空间W=[w1,w2,,wm],那么可以得到xi在空间W上的坐标为:Y=WTx∈ℜm*1。 xθ证明:wTxi=w⋅x⋅cosθ wy 又∵x⋅cosθ=y, w=1 ⇒wTxi=y 图 1 投影图 进一步,表达式w=m+ae表示w是一条通过点m,方向为e的直线。 1.2.2 PCA的作用及其统计特性[3] 采用PCA对原始数据的处理,通常有三个方面的作用—降维、相关性去除、概率估计。下面分别进行介绍: 去除原始数据相关性 从统计学上讲,E{[X−E(X)][Y−E(Y)]}称为随机变量X与Y协方差,记为 Cov(X,Y)。令ρXY= Cov(X,Y)D(X)D(Y),称为随机变量X与Y的相关系数。ρXY=1则X与 Y是相关的,ρXY=0,则X与Y是不相关的。 命题 1对于矩阵A来说,如果AA是一个对角阵,那么A中的向量是非相关的。 T 由PCA处理的人脸库数据的非相关性可以从两点进行说明。 (1) 基底的非相关性 特征空间基U=[u0,u1,,ur−1]是非相关的,即UUT=I。 (2) 投影系数的非相关性 根由SVD可知A={φ1,φ2,...,φM}=UΛVT, 其中φi=xi−mx,mx是平均人脸。据公式(2)可以把A映射到特征空间上,得到: B=UT*A,其中B是非相 关的,可由下面得到证明: Y的协方差矩阵为:CB= 1112 BBT=UTAATU=Λ (3) MMM 由命题 1可知,B是非相关的。 统计参数(均值及方差) 均值即mx--平均人脸。 命题 2随机变量方差越大,包含的信息越多,当一个变量方差为0时,该变量为常数,不含任何信息。 用PCA计算主分量,就是寻找一组向量,使得原始数据A={φ1,φ2,...,φM}在这组向量上的投影值的方差尽可能大。最大方差对应的向量就是第一主成份,以后递推就是第二主成份,第三主成份……。 用PCA计算主分量就是求原始数据A={φ1,φ2,...,φM}(其中φi=xi−mx)协方差矩阵的特征向量U=[u0,u1,,ur−1],由公式(3)可知,P=uiTA=(p1,p2,,pm)是A在ui上的投影值,其中P的方差就是ui对应的特征值λi,可以理解为: 命题 3 所有原始数据在主分量ui上的投影值方差为λi。 降维 如果在原始空间表示一幅n*m大小的图片X,那么需要一个N=n*m维矢量,但是当用公式(2)把它映射到特征空间后,只需要一个r*1维的向量就可。 另外,由命题2可知,可以根据方差的大小来判断特征向量的重要性。由ORL图片库的200个人脸计算得到的特征值呈图 2分布,可知特征向量重要性呈指数下降,据此可以只选用前面几个重要的特征向量来构建特征空间。 通过计算,前71个特征值占了90.17%,因此r可以取71而非200,从而达到进一步降维的作用。 图 2 特征值的分布 1.2.3 特征脸 U=[u0,u1,,ur−1]中的每一个单位向量都构成一个特征脸,如图 3所示。由这些特征脸 所张成的空间称为特征脸子空间,需要注意对于正交基的选择的不同考虑,对应较大特征值的特征向量(正交基)也称主分量,用于表示人脸的大体形状,而对应于较小特征值的特征向量则用于描述人脸的具体细节,或者从频域来看, 主分量表示了人脸的低频部分,而此分量则描述了人脸的高频部分(源码见’EigenFace.m’)。 1 2 10 50 70 average 图 3 特征脸,分别是第1,2,10,50,70分量,最后一张是平均脸。 1.2.4 图片重建 要进行图片X的重建,首先对X投影到特征空间上,得到系数Y=UT(X−mX),然后选用一部分系数与特征向量进行原始图片的重建:X'=mX+U(1:t)*Y(1:t),其中1:t表示取前t个元素。(见’reconstruct.m’) 在图 4中,其中前两张图片来自训练样本,第3张来自测试样本,可以看到对于训练样本,PCA系数可以对图片实现很好重建,而对于训练样本以外的图片重建效果很差。 Original1550100150199 图 4 人脸图像重建。第列张图片是输入原始图,其它列图片是重建结果,数字表示t的数目。 1.2.5 奇异值分解(SVD)[1] 设A是秩为r的m*n(m>>n)维矩阵,则存在两个正交矩阵和一个对角阵: A=[a1,a2,,ar]=UΛVT 其中U=[u0,u1,,ur−1],V=[v0,v1,,vr−1],Λ=diag(λ0,λ1,,λr−1),且UUT=I,VVT=I, λi呈降序排列。其中λi2为AAT∈ℜm*m和ATA∈ℜn*n的非零特征值,ui和vi分别是AAT和ATA对应于λi2的特征向量。可得一个推论: U=AVΛ−1 可以计算ATA的特征值λi2及相应的正交归一特征向量vi后,可由推论知AAT的正交归一特征向量 ui= 1 λi Avi 注意,协方差矩阵CA=(AiAT)/M的特征值为:λi2/M。 1.2.6 利用小矩阵计算大矩阵特征向量 高阶矩阵的特征向量可以转化为求低阶矩阵的特征向量: 设:A是秩为r 的m*n(m>>n)维矩阵,CX=AAT∈ℜm*m,是一个矩阵,现在要求CX的特征值及特征向量,可通过先求小矩阵ATA∈ℜn*n的特征向量[v0,v1,,vr−1]和特征值 [λ0,λ1,,λr−1],两者之间有以下关系: ATA⋅vi=λi⋅vi 左乘A ⎯⎯⎯→AAT(A⋅vi)=λi(A⋅vi) 显然,CX=AAT的特征向量是A⋅vi(注意没有单位化),[λ0,λ1,,λr−1]亦为其特征值。 结论:1.2.5与1.2.6的方法计算协方差矩阵的特征向量,特征值的结果是一致的,只是要注意1.2.5中的特征值要除以M,1.2.6中的特征向量要单位化。 1.2.7 图片归一化 图片标准化通常是一个整体概念,要求把图片归一到均值为0,方差为1下情况下。 这个概念类似于一般正态分布向标准正态分布的转化: 命题 4 若X∼N(µ,σ2),则Z= X−µσ所以要对一组图片中的一张Xi进行归一化(标准化),只需要减去均值,除以方差就可以了。 T ⎤均值mX=∑XiM,方差为D=E⎡(Xm)(Xm)−−XX⎣⎦ iM ∼N(0,1) 1.3 参考文献 [1] 邓楠, 基于主成份分析的人脸识别, 西北大学硕士学位论文,2006.06 [2] R.O. Duda, P.E. Hart, and D.G. Stork, Pattern Classification, seconded. John Wiley & Sons, 2001. [3] Sami Romdhani .Face Image Analysis using a Multiple Feature Fitting Strategy. PhD Thesis, University of Basel, Switzerland, January 2005. [4] Sami Romdhani .FACE RECOGNITION USING PRINCIPAL COMPONENTS ANALYSIS. 1.4 附录—matlab源码 1.4.1 人脸识别 % FaceRec.m % PCA 人脸识别修订版,识别率88% % calc xmean,sigma and its eigen decomposition allsamples=[];%所有训练图像 for i=1:40 for j=1:5 a=imread(strcat('e:\\ORL\\s',num2str(i),'\\',num2str(j),'.jpg')); % imshow(a); 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 for 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)); end i=1; % (训练阶段)计算特征脸形成的坐标系 base = xmean' * vsort(:,1:p) * diag(dsort(1:p).^(-1/2)); % base是N×p阶矩阵,除以dsort(i)^(1/2)是对人脸图像的标准化(使其方差为1) % 详见《基于PCA的人脸识别算法研究》p31 % 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) % 详见《基于PCA的人脸识别算法研究》p31 % i = i + 1; % xmean' * vsort(:,i)是小矩阵的特征向量向大矩阵特征向量转换的过程 %end % 以下两行add by gongxun 将训练样本对坐标系上进行投影,得到一个 M*p 阶矩阵allcoor allcoor = allsamples * base; % allcoor 里面是每张训练人脸图片在M*p子空间中的一个点,即在子空间中的组合系数, accu = 0; % 下面的人脸识别过程中就是利用这些组合系数来进行识别 % 测试过程 for i=1:40 for j=6:10 %读入40 x 5 副测试图像 a=imread(strcat('e:\\ORL\\s',num2str(i),'\\',num2str(j),'.jpg')); 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)-1)/5 )+1; class2=floor((index2(2)-1)/5)+1; class3=floor((index2(3)-1)/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 %输出识别率 1.4.2 特征人脸 % eigface.m function [] = eigface() % calc xmean,sigma and its eigen decomposition allsamples=[];%所有训练图像 for i=1:40 for j=1:5 a=imread(strcat('e:\\ORL\\s',num2str(i),'\\',num2str(j),'.jpg')); % imshow(a); 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 for 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)); end p = 199; % (训练阶段)计算特征脸形成的坐标系 %while (i<=p && dsort(i)>0) % base(:,i) = dsort(i)^(-1/2) * xmean' * vsort(:,i); % base是N×p阶矩阵,除以dsort(i)^(1/2)是对人脸图像的标准化,详见《基于PCA的人脸识别算法研究》p31 % i = i + 1; % xmean' * vsort(:,i)是小矩阵的特征向量向大矩阵特征向量转换的过程 %end base = xmean' * vsort(:,1:p) * diag(dsort(1:p).^(-1/2)); % 生成特征脸 for (k=1:p), temp = reshape(base(:,k), 112,92); newpath = ['d:\est\\' int2str(k) '.jpg']; imwrite(mat2gray(temp), newpath); end avg = reshape(samplemean, 112,92); imwrite(mat2gray(avg), 'd:\est\\average.jpg'); % 将模型保存 save('e:\\ORL\\model.mat', 'base', 'samplemean'); 1.4.3 人脸重建 % Reconstruct.m function [] = reconstruct() load e:\\ORL\\model.mat; % 计算新图片在特征子空间中的系数 img = 'D:\est2\\10.jpg' a=imread(img); b=a(1:112*92); % b是行矢量 1×N,其中N=10304,提取顺序是先列后行,即从上到下,从左到右 b=double(b); b=b-samplemean; c = b * base; % c 是图片a在子空间中的系数, 是1*p行矢量 % 根据特征系数及特征脸重建图 % 前15个 t = 15; temp = base(:,1:t) * c(1:t)'; temp = temp + samplemean'; imwrite(mat2gray(reshape(temp, 112,92)),'d:\est2\1.jpg'); % 前50个 t = 50; temp = base(:,1:t) * c(1:t)'; temp = temp + samplemean'; imwrite(mat2gray(reshape(temp, 112,92)),'d:\est2\2.jpg'); % 前100个 t = 100; temp = base(:,1:t) * c(1:t)'; temp = temp + samplemean'; imwrite(mat2gray(reshape(temp, 112,92)),'d:\est2\3.jpg'); % 前150个 t = 150; temp = base(:,1:t) * c(1:t)'; temp = temp + samplemean'; imwrite(mat2gray(reshape(temp, 112,92)),'d:\est2\4.jpg'); % 前199个 t = 199; temp = base(:,1:t) * c(1:t)'; temp = temp + samplemean'; imwrite(mat2gray(reshape(temp, 112,92)),'d:\est2\5.jpg'); 因篇幅问题不能全部显示,请点此查看更多更全内容