EigenFace算法詳解及Matlab代碼2014

1、EigenFace算法回顧及Matlab代碼 基于PCA的EigenFace算法發(fā)表自1987年，是第一種可行的人臉辨識(shí)算法。雖然已有20余年歷史，但仍是人臉辨識(shí)算法研究中的經(jīng)典，新算法都要與之作比較。EigenFace是2D辨識(shí)算法，但為了進(jìn)行3D表情辨識(shí)的研究，有必要對(duì)這一經(jīng)典算法進(jìn)行回顧，SIGGRAPH13的文獻(xiàn)OnlineModelingForRealtimeFacialAnimation實(shí)現(xiàn)表情3D重構(gòu)的基礎(chǔ)是SIGGRAPH99中AMorphableModelfortheSynthesisof3DFaces提出的MorphableFacialModel，而建立這一模型的基礎(chǔ)思想仍

2、是PCA，與Eigenface有著天然聯(lián)系。學(xué)習(xí)EigenFace應(yīng)該是研究生時(shí)代的事兒了，舊編重拾、開(kāi)卷有益，并寫了Matlab代碼附錄于后。網(wǎng)上許多實(shí)例代碼只實(shí)現(xiàn)了辨識(shí)，略去了一個(gè)重要環(huán)節(jié):通過(guò)分解-重構(gòu)，將一幅輸入人臉照片表示為EigenFace基底的組合，這對(duì)于表情辨識(shí)及3D人臉模型分析都是很有用的(OnlineModelingForRealtimeFacialAnimation中就利用了此思想)，附錄的Matlab代碼做了這一步。訓(xùn)練集合包含20幅圖片生成20-1=19個(gè)特征臉 最近歐式距離法得出的前三位匹配 利用特征臉空間進(jìn)行人臉重構(gòu)http:/www.pages.drexel.e

3、du/sis26/EigenfaceTutorial.htm中給出了EigenFace算法Matlab代碼，含有重構(gòu)過(guò)程與比較，但其代碼中有一個(gè)錯(cuò)誤，訓(xùn)練集合照片未減去平均臉，就計(jì)算協(xié)方差矩陣了具體算法文獻(xiàn)教材上都有，撿要點(diǎn)寫幾句：(1)將訓(xùn)練集合中的每幅圖像拉伸為列向量，并減去所有圖像的均值(稱之為平均臉)，形成N*M矩陣A,其中N為單幅圖像像素?cái)?shù)，M為圖像數(shù)目(訓(xùn)練集合容量)。(2)求協(xié)方差矩陣AA的特征向量，作為正交基底張成人臉空間，好是好但運(yùn)算量過(guò)大，轉(zhuǎn)而求替代矩陣(surrogate)AA的特征向量，減少計(jì)算量(3)矩陣AA的秩等于M-1，這是由于減去平均臉?biāo)?，故有M-1個(gè)非零特征

4、值(正)，去除屬于0的特征向量，將M-1個(gè)屬于非零特征值得特征向量(記住須作左乘A的修正)作為EigenFace基底(特征臉)，張成人臉空間。(4)EigenFace基底由M-1個(gè)相互正交的向量構(gòu)成，它們是協(xié)方差矩陣AA的前M-1個(gè)最顯著的特征向量方向，能量主要集中在這些向量方向上，但要記住雖正交但不完備，故存在重構(gòu)誤差。(5)為了用EigenFace基底對(duì)人臉照片進(jìn)行正確的分解-重構(gòu)，需要對(duì)所得的基底向量進(jìn)行規(guī)一化修正，因?yàn)锳A的特征向量左乘A之后，雖成為AA的特征向量，但模不為1，需除以自身的模，修正為標(biāo)準(zhǔn)正交向量集合，才能進(jìn)行投影分解-重構(gòu)。(6)訓(xùn)練集合及測(cè)試集合中的人臉照片都能利用E

5、igenFace基底較好地實(shí)現(xiàn)分解-重構(gòu)，但訓(xùn)練集合之外的人臉，重構(gòu)誤差變大(7)EigenFace缺點(diǎn)：1拍攝時(shí)光照環(huán)境對(duì)識(shí)別效果(EigenFace基底)影響大2訓(xùn)練集合擴(kuò)容時(shí)，需重構(gòu)EigenFace基底(8)為了有效顯示EigenFace基底圖像(特征臉)需要用imagesc函數(shù)Eigenface算法識(shí)別人臉的步驟：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 . 新臉向量在每個(gè)特征臉向量上的投影：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代碼：所用訓(xùn)練集合含20幅jpg圖片，存于名為TrainDataBase的文件夾中，測(cè)試集合含10幅jpg圖片，存于名為TestDataBase的文件夾中，命名均為1.jpg2.jpg3.jpg.%-%人臉辨識(shí)：Face Recognition%EigenFace算法%EigenFace空間的生成

16、分解與合成% 人臉辨識(shí)%-function meigenface()%-%step1:讀取訓(xùn)練集合中的人臉照片，并通過(guò)列向量拉伸生成% 人臉訓(xùn)練集合矩陣%-clear allclcclose allimset=;M=20;%訓(xùn)練集合中人臉照片的數(shù)目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中數(shù)據(jù)轉(zhuǎn)化為imheight*imwidth高的列向量 imset=imset imv; % 將每個(gè)圖像的像素占據(jù)一列位置，共20列向量if i=4title(訓(xùn)練集合人臉照片,fontsize,18)endsubplot(5,4,i);imshow(im);enddisp(=Output image is completed!=)% imset=double(imset);imset=double(imset);w,h=size(imset)%-%step2:生成訓(xùn)練集合的平均臉%-averageface=mean(imset,2);% imset每行的平均值，

18、由于imset為單列向量，每行包含一個(gè)元素，因此此操作為將每行的所有值取平均值figure(2);imshow(uint8(reshape(averageface,imheight,imwidth);title(平均臉);%-%step3:將人臉訓(xùn)練集合修正為與平均人臉的差:以便計(jì)算% 協(xié)方差矩陣 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特征臉空間%-%互相關(guān)矩陣為A*A,為減少計(jì)算量，計(jì)算代理矩陣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(zhǎng)的特征向量，d為L(zhǎng)的特征值dL=L*v-v*d%人臉照片修正為與平均臉的差，此時(shí)矩陣A*A的秩為M-1,特征值中含有一個(gè)0%需要去除屬于0的特征向量,Eigenface基底由M-1組正交向量組成v d= EliminateZero(v,d) %去除屬于0的特征向量，同時(shí)將特征向量按特征值大小降序排列%生成人臉空間矩陣:eigenfaceeigenface=;for i=1:M-1mv=A*v(:,i);%左乘矩陣A，從而修正為協(xié)方差矩陣AA的特征向量%對(duì)Eigenface向量進(jìn)行規(guī)一化修正，使每個(gè)Eigenface向量模為1，便

21、于進(jìn)行分解與合成%是否進(jìn)行此歸一化修正對(duì)識(shí)別沒(méi)有影響，但對(duì)利用Eigenface作為基底對(duì)人臉照片進(jìn)行%分解與重構(gòu)卻是必要的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: 計(jì)算訓(xùn)練集合中的照片在EigenFace基底上的投影坐標(biāo)% 以便進(jìn)行人臉辨識(shí)%注:如果用Eigenface基底直接對(duì)訓(xùn)練集合中的照片進(jìn)行分解-重構(gòu)%則重構(gòu)誤差趨于0(1e-10)，如果對(duì)測(cè)試集合中的人臉(同一人，但表情姿態(tài)略有出入)%進(jìn)行重構(gòu)，則可觀察到重構(gòu)誤差，這是由于EigenFace基底雖正交但不完備%- TrainCoeff=;for i=1:M%計(jì)算訓(xùn)練集合中第i幅圖像在EigenFace基底上的投影disp(The %d th face image)i=ic=do

23、uble(imset_d(:,i)*eigenface;TrainCoeff=TrainCoeff cend%-%step6:人臉辨識(shí)，采用最近歐式距離法%計(jì)算輸入圖像在EigenFace基底上的投影坐標(biāo)%與訓(xùn)練集合中各圖像的投影坐標(biāo)進(jìn)行比較，選取歐式距離最小者作為匹配項(xiàng)%- %讀入測(cè)試集合人臉照片N=10;%測(cè)試集合人臉照片數(shù)目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,請(qǐng)選擇一幅照片(1-10),fontsize,18);testset=double(testset);n= input(請(qǐng)選擇一幅照片(1-10): ); %選取訓(xùn)練集合中第n幅人臉照片作為輸入進(jìn)行分解InputImage=testset(:,n)-averageface;%計(jì)算輸入照片在EigenFace基底上的投影坐標(biāo)向量coeff=Input

25、Image*eigenface;%計(jì)算與訓(xùn)練集合中各圖像的投影坐標(biāo)向量的歐式距離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基底對(duì)測(cè)試集合中的人臉照片進(jìn)行% 分解-合成(重構(gòu))%注:如果用Eigenface基底直接對(duì)訓(xùn)練集合中的照片進(jìn)行分解-重構(gòu)%則重構(gòu)誤差趨于0(1e-10)，如果對(duì)測(cè)試集合中的人臉(同一人，但表情姿態(tài)略有出入)%進(jìn)行重構(gòu)，則可觀察到重構(gòu)誤差，這是由于EigenFace基底雖正交但不完備%- ;%選取訓(xùn)練集合中第n幅人臉照片作為輸入進(jìn)行分解InputImage=testset(:,n)-averageface;%計(jì)算輸入照片在EigenFace基底上的投影坐標(biāo)coeff=InputImage*eigenface;%在EigenFace基底上進(jìn)行圖像重構(gòu)Recon

28、structedImage=eigenface(:,1:M-1)*coeff;%計(jì)算重構(gòu)誤差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(訓(xùn)練集合中照片1);subplot(2,2,2);im2=imset(:,n*2);im2=reshape(im2,imheight,imwidth);imshow(uint8(im2);title(訓(xùn)練集合中照片2);subplot(2,2,3);InputImage=reshape(InputImage,imheight,imwidth);imshow(uint8(InputImage);title(測(cè)試照片);subp

30、lot(2,2,4);ReconstructedImage=reshape(ReconstructedImage,imheight,imwidth);imshow(uint8(ReconstructedImage);title(測(cè)試照片的重構(gòu));pause;return;%-%去除屬于0的特征向量，同時(shí)將特征向量按特征值大小降序排列%-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=