EigenFace算法详解及Matlab代码2014

1、EigenFace算法回顾及Matlab代码 基于PCA的EigenFace算法发表自1987年，是第一种可行的人脸辨识算法。虽然已有20余年历史，但仍是人脸辨识算法研究中的经典，新算法都要与之作比较。EigenFace是2D辨识算法，但为了进行3D表情辨识的研究，有必要对这一经典算法进行回顾，SIGGRAPH13的文献OnlineModelingForRealtimeFacialAnimation实现表情3D重构的基础是SIGGRAPH99中AMorphableModelfortheSynthesisof3DFaces提出的MorphableFacialModel，而建立这一模型的基础思想仍

2、是PCA，与Eigenface有着天然联系。学习EigenFace应该是研究生时代的事儿了，旧编重拾、开卷有益，并写了Matlab代码附录于后。网上许多实例代码只实现了辨识，略去了一个重要环节:通过分解-重构，将一幅输入人脸照片表示为EigenFace基底的组合，这对于表情辨识及3D人脸模型分析都是很有用的(OnlineModelingForRealtimeFacialAnimation中就利用了此思想)，附录的Matlab代码做了这一步。训练集合包含20幅图片生成20-1=19个特征脸 最近欧式距离法得出的前三位匹配 利用特征脸空间进行人脸重构http:/www.pages.drexel.e

3、du/sis26/EigenfaceTutorial.htm中给出了EigenFace算法Matlab代码，含有重构过程与比较，但其代码中有一个错误，训练集合照片未减去平均脸，就计算协方差矩阵了具体算法文献教材上都有，捡要点写几句：(1)将训练集合中的每幅图像拉伸为列向量，并减去所有图像的均值(称之为平均脸)，形成N*M矩阵A,其中N为单幅图像像素数，M为图像数目(训练集合容量)。(2)求协方差矩阵AA的特征向量，作为正交基底张成人脸空间，好是好但运算量过大，转而求替代矩阵(surrogate)AA的特征向量，减少计算量(3)矩阵AA的秩等于M-1，这是由于减去平均脸所致，故有M-1个非零特征

4、值(正)，去除属于0的特征向量，将M-1个属于非零特征值得特征向量(记住须作左乘A的修正)作为EigenFace基底(特征脸)，张成人脸空间。(4)EigenFace基底由M-1个相互正交的向量构成，它们是协方差矩阵AA的前M-1个最显著的特征向量方向，能量主要集中在这些向量方向上，但要记住虽正交但不完备，故存在重构误差。(5)为了用EigenFace基底对人脸照片进行正确的分解-重构，需要对所得的基底向量进行规一化修正，因为AA的特征向量左乘A之后，虽成为AA的特征向量，但模不为1，需除以自身的模，修正为标准正交向量集合，才能进行投影分解-重构。(6)训练集合及测试集合中的人脸照片都能利用E

5、igenFace基底较好地实现分解-重构，但训练集合之外的人脸，重构误差变大(7)EigenFace缺点：1拍摄时光照环境对识别效果(EigenFace基底)影响大2训练集合扩容时，需重构EigenFace基底(8)为了有效显示EigenFace基底图像(特征脸)需要用imagesc函数Eigenface算法识别人脸的步骤：This section gives step-by-step instructions along with photos and formulas on how to recognize faces and implemented into Matlab. All th

6、e necessary files to complete this tutorial would be provided.Steps1 The first step is to obtain a set S with M face images. In our example M = 25 as shown at the beginning of the tutorial. Each image is transformed into a vector of size N and placed into the set.2 After you have obtained your set,

7、you will obtain the mean image 3 Then you will find the difference between the input image and the mean image4 Next we seek a set of M orthonormal vectors, un, which best describes the distribution of the data. The kth vector, uk, is chosen such that 求特征值：is a maximum, subject toNote: uk and k are t

8、he eigenvectors and eigenvalues of the covariance matrix C 5 We obtain the covariance matrix C in the following manner Where the matrix . The matrix C, however, is N2 by N2, and determining the N2 eigenvectors and eigenvalues is an intractable task for typical image sizes. We need a computationally

9、feasible method to find these eigenvectors. We can solve for the N2 dimensional eigenvectors in this case by first solving for the eigenvectors of an M by M matrix- e.g., solving a 16*16 matrix rather than a 16,384*16,384 matrix and then taking appropriate linear combinations of the face images i,Co

10、nsider the eigenvectors vi of ATA such that ATAvi=iviPremultiplying both sides by A, we have AATAvi=iAviFrom which we see that Avi (特征向量)are the eigenvectors of C=AAT.Following this analysis, we construct the M by M matrix L=ATA, where(是1*16,384维的向量), and find the M eigenvectors, vi (M维的向量)of L. The

11、se vectors determine linear combinations of the M training set face images to form the eigenfaces ul (特征脸向量). Once we have found the eigenvectors, vl, ul ul=1,2,Mvlt=Avlt=k=1Mvlkk l=1,2,M6 These are the eigenfaces of our set of original images Recognition Procedure1. A new face is transformed into i

12、ts eigenface components. First we compare our input image with our mean image and multiply their difference with each eigenvector of the L matrix. Each value would represent a weight and would be saved on a vector . 新脸向量在每个特征脸向量上的投影：For k=1,2,M. 2. We now determine which face class provides the best

13、 description for the input image. This is done by minimizing the Euclidean distance 3. The input face is consider to belong to a class if k is bellow an established threshold . Then the face image is considered to be a known face. If the difference is above the given threshold, but bellow a second t

14、hreshold, the image can be determined as a unknown face. If the input image is above these two thresholds, the image is determined NOT to be a face. 4. If the image is found to be an unknown face, you could decide whether or not you want to add the image to your training set for future recognitions.

15、 You would have to repeat steps 1 trough 7 to incorporate this new face image. The source code face recognition using Matlab is provided below: Matlab代码：所用训练集合含20幅jpg图片，存于名为TrainDataBase的文件夹中，测试集合含10幅jpg图片，存于名为TestDataBase的文件夹中，命名均为1.jpg2.jpg3.jpg.%-%人脸辨识：Face Recognition%EigenFace算法%EigenFace空间的生成

16、分解与合成% 人脸辨识%-function meigenface()%-%step1:读取训练集合中的人脸照片，并通过列向量拉伸生成% 人脸训练集合矩阵%-clear allclcclose allimset=;M=20;%训练集合中人脸照片的数目figure(1);for i=1:Mfilename=.TrainDataBase num2str(i) .jpg;im=imread(filename,JPG);im=im(:,:,2);% im=im(:,:,3);imheight imwidth=size(im);%将图像矩阵按列拉伸为列向量imv=reshape(im,imheight*i

17、mwidth,1); %将im中数据转化为imheight*imwidth高的列向量 imset=imset imv; % 将每个图像的像素占据一列位置，共20列向量if i=4title(训练集合人脸照片,fontsize,18)endsubplot(5,4,i);imshow(im);enddisp(=Output image is completed!=)% imset=double(imset);imset=double(imset);w,h=size(imset)%-%step2:生成训练集合的平均脸%-averageface=mean(imset,2);% imset每行的平均值，

18、由于imset为单列向量，每行包含一个元素，因此此操作为将每行的所有值取平均值figure(2);imshow(uint8(reshape(averageface,imheight,imwidth);title(平均脸);%-%step3:将人脸训练集合修正为与平均人脸的差:以便计算% 协方差矩阵 Covariance matrix%-imset_d=imset;figure(3);for i=1:20imset_d(:,i)=imset(:,i)- averageface;im=reshape(imset_d(:,i),imheight,imwidth);if i=4title(与平均脸的差

19、,fontsize,18)end% subplot(4,5,i);subplot(5,4,i);im=imagesc(im); colormap(gray);% imshow(im);end% %-%step4:生成EigenFace特征脸空间%-%互相关矩阵为A*A,为减少计算量，计算代理矩阵A*A的特征值特征向量A=double(imset_d); % include the data of twenty images. % matrix of 36000 hight and 20 widthAh,Aw=size(A)L=A*A; % matrix of 20 hight and 20 w

20、idthLh,Lw=size(L)v d=eig(L) % v为L的特征向量，d为L的特征值dL=L*v-v*d%人脸照片修正为与平均脸的差，此时矩阵A*A的秩为M-1,特征值中含有一个0%需要去除属于0的特征向量,Eigenface基底由M-1组正交向量组成v d= EliminateZero(v,d) %去除属于0的特征向量，同时将特征向量按特征值大小降序排列%生成人脸空间矩阵:eigenfaceeigenface=;for i=1:M-1mv=A*v(:,i);%左乘矩阵A，从而修正为协方差矩阵AA的特征向量%对Eigenface向量进行规一化修正，使每个Eigenface向量模为1，便

21、于进行分解与合成%是否进行此归一化修正对识别没有影响，但对利用Eigenface作为基底对人脸照片进行%分解与重构却是必要的mv=mv/norm(mv); %!eigenface=eigenface mv; % 特征人脸空间eh,ew=size(eigenface)end%显示归一化之后的人脸空间eigenfacefigure(4);for i=1:M-1im=eigenface(:,i);im=reshape(im,imheight,imwidth);if i=4title(特征脸空间：EigenFace,fontsize,18)endsubplot(5,4,i);im=imagesc(im

22、);colormap(gray);%im=histeq(im,255);imshow(im);end%-%step5: 计算训练集合中的照片在EigenFace基底上的投影坐标% 以便进行人脸辨识%注:如果用Eigenface基底直接对训练集合中的照片进行分解-重构%则重构误差趋于0(1e-10)，如果对测试集合中的人脸(同一人，但表情姿态略有出入)%进行重构，则可观察到重构误差，这是由于EigenFace基底虽正交但不完备%- TrainCoeff=;for i=1:M%计算训练集合中第i幅图像在EigenFace基底上的投影disp(The %d th face image)i=ic=do

23、uble(imset_d(:,i)*eigenface;TrainCoeff=TrainCoeff cend%-%step6:人脸辨识，采用最近欧式距离法%计算输入图像在EigenFace基底上的投影坐标%与训练集合中各图像的投影坐标进行比较，选取欧式距离最小者作为匹配项%- %读入测试集合人脸照片N=10;%测试集合人脸照片数目figure(5);testset=;for i=1:Nfilename=.TestDataBase num2str(i) .jpg;im=imread(filename,JPG);im=im(:,:,2);imheight imwidth=size(im);%将图像

24、矩阵按列拉伸为列向量imv=reshape(im,imheight*imwidth,1);testset=testset imv;% subplot(3,4,i);imshow(im);title(int2str(i);endtext(200,15,请选择一幅照片(1-10),fontsize,18);testset=double(testset);n= input(请选择一幅照片(1-10): ); %选取训练集合中第n幅人脸照片作为输入进行分解InputImage=testset(:,n)-averageface;%计算输入照片在EigenFace基底上的投影坐标向量coeff=Input

25、Image*eigenface;%计算与训练集合中各图像的投影坐标向量的欧式距离cdist=;for i=1:Mdist=(norm(coeff-TrainCoeff(:,i)2;cdist=cdist dist;endd index=sort(cdist)inum=index(1),index(2),index(3);figure(6);subplot(2,3,2);InputImage=reshape(InputImage+averageface,imheight,imwidth);imshow(uint8(InputImage);title(输入照片);subplot(2,3,4);im

26、1=imset(:,index(1);im1=reshape(im1,imheight,imwidth);imshow(uint8(im1);title(匹配照片1);subplot(2,3,5);im2=imset(:,index(2);im2=reshape(im2,imheight,imwidth);imshow(uint8(im2);title(匹配照片2);subplot(2,3,6);im2=imset(:,index(3);im2=reshape(im2,imheight,imwidth);imshow(uint8(im2);title(匹配照片3);pause;close;%-

27、%step6:利用EigenFace基底对测试集合中的人脸照片进行% 分解-合成(重构)%注:如果用Eigenface基底直接对训练集合中的照片进行分解-重构%则重构误差趋于0(1e-10)，如果对测试集合中的人脸(同一人，但表情姿态略有出入)%进行重构，则可观察到重构误差，这是由于EigenFace基底虽正交但不完备%- ;%选取训练集合中第n幅人脸照片作为输入进行分解InputImage=testset(:,n)-averageface;%计算输入照片在EigenFace基底上的投影坐标coeff=InputImage*eigenface;%在EigenFace基底上进行图像重构Recon

28、structedImage=eigenface(:,1:M-1)*coeff;%计算重构误差diff=ReconstructedImage-InputImage;norm(diff),InputImage=reshape(InputImage+averageface,imheight,imwidth);ReconstructedImage=reshape(ReconstructedImage+averageface,imheight,imwidth);%figure(5);%imagesc(ReconstructedImage); colormap(gray);figure(7);subplo

29、t(2,2,1);im1=imset(:,n*2-1);im1=reshape(im1,imheight,imwidth);imshow(uint8(im1);title(训练集合中照片1);subplot(2,2,2);im2=imset(:,n*2);im2=reshape(im2,imheight,imwidth);imshow(uint8(im2);title(训练集合中照片2);subplot(2,2,3);InputImage=reshape(InputImage,imheight,imwidth);imshow(uint8(InputImage);title(测试照片);subp

30、lot(2,2,4);ReconstructedImage=reshape(ReconstructedImage,imheight,imwidth);imshow(uint8(ReconstructedImage);title(测试照片的重构);pause;return;%-%去除属于0的特征向量，同时将特征向量按特征值大小降序排列%-function v d= EliminateZero(v,d)%取特征值% dd=diag(d);dd=abs(diag(d);%找出等于0的特征值ind=find(dd1e-06);% 1e-06);% /1E-06);% 1e-06);% 1e-06);% /1E-06);%去除=0的特征值dd(ind)=;%去除等于0的特征向量v(:,ind)=;%特征值降序排序ddd ind=sort(dd,descend);d=