【老生谈算法】matlab实现EigenFace算法源码——EigenFace算法

EigenFace算法回顾及Matlab代码

1、文档下载：

本算法已经整理成文档如下，有需要的朋友可以点击进行下载

2、算法详解：

基于PCA的EigenFace算法发表自1987年，是第一种可行的人脸辨识算法。虽然已有20余年历史，但仍是人脸辨识算法研究中的经典，新算法都要与之作比较。

EigenFace是2D辨识算法，但为了进行3D表情辨识的研究，有必要对这一经典算法进行回顾，SIGGRAPH13的文献Online Modeling For Realtime Facial Animation实现表情3D重构的基础是SIGGRAPH99中A Morphable Model for the Synthesis of 3D Faces提出的Morphable Facial Model，而建立这一模型的基础思想仍是PCA，与Eigenface有着天然联系。

学习EigenFace应该是研究生时代的事儿了，旧编重拾、开卷有益，并写了Matlab代码附录于后。网上许多实例代码只实现了辨识，略去了一个重要环节:通过分解-重构，将一幅输入人脸照片表示为EigenFace基底的组合，这对于表情辨识及3D人脸模型分析都是很有用的(Online Modeling For Realtime Facial Animation中就利用了此思想)，附录的Matlab代码做了这一步。

                                      训练集合包含20幅图片

生成20-1=19个特征脸

最近欧式距离法得出的前三位匹配

利用特征脸空间进行人脸重构
http://www.pages.drexel.edu/~sis26/Eigenface Tutorial.htm中给出了EigenFace算法Matlab 代码，含有重构过程与比较，但其代码中有一个错误，训练集合照片未减去平均脸，就计算协方差矩阵了

具体算法文献教材上都有，捡要点写几句：
(1) 将训练集合中的每幅图像拉伸为列向量，并减去所有图像的均值(称之为平均脸)，形成N*M矩阵A,其中N为单幅图像像素数，M为图像数目(训练集合容量)。
(2) 求协方差矩阵AA’的特征向量，作为正交基底张成人脸空间，好是好但运算量过大，转而求替代矩阵(surrogate)A’A的特征向量，减少计算量
(3) 矩阵A’A的秩等于M-1，这是由于减去平均脸所致，故有M-1个非零特征值(正)，去除属于0的特征向量，将M-1个属于非零特征值得特征向量(记住须作左乘A的修正)作为EigenFace基底(特征脸)，张成人脸空间。
(4) EigenFace基底由M-1个相互正交的向量构成，它们是协方差矩阵AA’的前M-1个最显著的特征向量方向，能量主要集中在这些向量方向上，但要记住虽正交但不完备，故存在重构误差。
(5) 为了用EigenFace基底对人脸照片进行正确的分解-重构，需要对所得的基底向量进行规一化修正，因为A’A的特征向量左乘A之后，虽成为AA’的特征向量，但模不为1，需除以自身的模，修正为标准正交向量集合，才能进行投影分解-重构。
(6) 训练集合及测试集合中的人脸照片都能利用EigenFace基底较好地实现分解-重构，但训练集合之外的人脸，重构误差变大
(7) EigenFace缺点：1拍摄时光照环境对识别效果(EigenFace基底)影响大
2训练集合扩容时，需重构EigenFace基底
(8) 为了有效显示EigenFace基底图像(特征脸)需要用imagesc函数

Matlab代码：
所用训练集合含20幅jpg图片，存于名为TrainDataBase的文件夹中，测试集合含10幅jpg图片，存于名为TestDataBase的文件夹中，命名均为1.jpg 2.jpg 3.jpg…

%----------------------------------------------------------------------
%人脸辨识：Face Recognition
%EigenFace算法
%EigenFace空间的生成分解与合成
%人脸辨识
%-----------------------------------------------------------------------
function meigenface()
%------------------------------------------------------
%step1:读取训练集合中的人脸照片，并通过列向量拉伸生成
%人脸训练集合矩阵
%------------------------------------------------------
clear all
clc
close all
imset=[];
M=20;%训练集合中人脸照片的数目
figure(1);
for i=1:M
filename=['.\TrainDataBase\' num2str(i) '.jpg'];
im=imread(filename,'JPG');
im=im(:,:,2);
% im=im(:,:,3);
[imheight imwidth]=size(im);
%将图像矩阵按列拉伸为列向量
imv=reshape(im,imheight*imwidth,1);%将im中数据转化为imheight*imwidth高的列向量
imset=[imset imv];%将每个图像的像素占据一列位置，共20列向量
if i==4
title('训练集合人脸照片','fontsize',18)
end
subplot(5,4,i);
imshow(im);
end
disp('============Output image is completed!===================')
% imset=double(imset);
imset=double(imset);
[w,h]=size(imset)
%------------------------------------------------------
%step2:生成训练集合的平均脸
%------------------------------------------------------
averageface=mean(imset,2);
% imset每行的平均值，由于imset为单列向量，每行包含一个元素，因此此操作为将每行的所有值取平均值
figure(2);
imshow(uint8(reshape(averageface,imheight,imwidth)));
title('平均脸');
%------------------------------------------------------
%step3:将人脸训练集合修正为与平均人脸的差:以便计算
%协方差矩阵Covariance matrix
%------------------------------------------------------
imset_d=imset;
figure(3);
for i=1:20
imset_d(:,i)=imset(:,i)- averageface;
im=reshape(imset_d(:,i),imheight,imwidth);
if i==4
title('与平均脸的差','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 width
[Ah,Aw]=size(A)
L=A'*A;   % matrix of 20 hight and 20 width
[Lh,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的特征向量，同时将特征向量按特征值大小降序排列
%生成人脸空间矩阵:eigenface
eigenface=[];
for i=1:M-1
mv=A*v(:,i);%左乘矩阵A，从而修正为协方差矩阵AA'的特征向量
%对Eigenface向量进行规一化修正，使每个Eigenface向量模为1，便于进行分解与合成
%是否进行此归一化修正对识别没有影响，但对利用Eigenface作为基底对人脸照片进行
%分解与重构却是必要的
mv=mv/norm(mv); %
eigenface=[eigenface mv];%特征人脸空间
[eh,ew]=size(eigenface)
end
%显示归一化之后的人脸空间eigenface
figure(4);
for i=1:M-1
im=eigenface(:,i);
im=reshape(im,imheight,imwidth);
if i==4
title('特征脸空间：EigenFace','fontsize',18)
end
subplot(5,4,i);
im=imagesc(im);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=i
c=double(imset_d(:,i)')*eigenface;
TrainCoeff=[TrainCoeff c']
end
%------------------------------------------------------
%step6:人脸辨识，采用最近欧式距离法
%计算输入图像在EigenFace基底上的投影坐标
%与训练集合中各图像的投影坐标进行比较，选取欧式距离最小者作为匹配项
%------------------------------------------------------
%读入测试集合人脸照片
N=10;%测试集合人脸照片数目
figure(5);
testset=[];
for i=1:N
filename=['.\TestDataBase\' num2str(i) '.jpg'];
im=imread(filename,'JPG');
im=im(:,:,2);
[imheight imwidth]=size(im);
%将图像矩阵按列拉伸为列向量
imv=reshape(im,imheight*imwidth,1);
testset=[testset imv];
%
subplot(3,4,i);
imshow(im);title(int2str(i));
end
text(200,15,'请选择一幅照片(1-10)','fontsize',18);
testset=double(testset);
n= input('请选择一幅照片(1-10): ');

%选取训练集合中第n幅人脸照片作为输入进行分解
InputImage=testset(:,n)-averageface;
%计算输入照片在EigenFace基底上的投影坐标向量
coeff=InputImage'*eigenface;
%计算与训练集合中各图像的投影坐标向量的欧式距离
cdist=[];
for i=1:M
dist=(norm(coeff'-TrainCoeff(:,i)))^2;
cdist=[cdist dist];
end
[d 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);
im1=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;
%------------------------------------------------------
%step6:利用EigenFace基底对测试集合中的人脸照片进行
%分解-合成(重构)
%注:如果用Eigenface基底直接对训练集合中的照片进行分解-重构
%则重构误差趋于0(1e-10)，如果对测试集合中的人脸(同一人，但表情姿态略有出入)
%进行重构，则可观察到重构误差，这是由于EigenFace基底虽正交但不完备
%------------------------------------------------------
;%选取训练集合中第n幅人脸照片作为输入进行分解
InputImage=testset(:,n)-averageface;
%计算输入照片在EigenFace基底上的投影坐标
coeff=InputImage'*eigenface;
%在EigenFace基底上进行图像重构
ReconstructedImage=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);
subplot(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('测试照片');
subplot(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(dd<1e-06);
% <1e-06);
% </1E-06);
% <1e-06);
% <1e-06);
% </1E-06);
%去除=0的特征值
dd(ind)=[];
%去除等于0的特征向量
v(:,ind)=[];
%特征值降序排序
[ddd ind]=sort(dd,'descend');
d=diag(ddd);
%重排特征向量
v=v(:,ind);
return;
/Result Pictures///







//
1.Example of reshape
tm=[]
ma=[1 9 7; 2 4 8;17 78 1; 22 4 79]
[imheight imwidth]=size(ma)
mac=reshape(ma,imheight*imwidth,1)
tm=[tm mac]
result:
mac =
1
2
17
22
9
4
78
4
7
8
1
79
2.Mean: compute mean value of matrix
mean(X,DIM) takes the mean along the dimension DIM of X.

Example: If X = [1 2 3; 3 3 6; 4 6 8; 4 7 7];

then mean(X,1) is [3.0000 4.5000 6.0000] and
mean(X,2) is [2.0000 4.0000 6.0000 6.0000].'
3.Imagesc
load clown
figure
subplot(1,2,1)
imshow(uint8(X))
subplot(1,2,2)
imagesc(X)
%colormap(gray)
axis image
title('Default CLim (= [1 81])')

4.Eig:Eigenvalues and eigenvectors.
E = eig(X) is a vector containing the eigenvalues of a square matrix X.
[V,D] = eig(X) produces a diagonal matrix D of eigenvalues and a
full matrix V whose columns are the corresponding eigenvectors so
that X*V = V*D.
B=[1 3 7; 31 8 4;11 89 3]
[VB,DB] = eig(B)
B*VB - VB*DB
[VN,DN] = eig(B,'nobalance')
B*VN - VN*DN
5.Norm
norm   Matrix or vector norm.
norm(X,2) returns the 2-norm of X.
norm(X) is the same as norm(X,2).