matlab实现pca降维并实现人脸识别

前言

题目复述

Experiment 7: PCA in Face Recognition

November 27, 2018

1 Description

The ORL Database of Faces (https://www.cl.cam.ac.uk/research/dtg/attarchive/ facedatabase.html) is a widely used database in the context of face recognition applications. The images are organized in 40 directories, each of which consists of 10 images. Note that, in each directories, the images are taken from the same person. To construct a set of training data, we randomly pick 5∼7 images from each categories, while leaving the others as test data. Please first apply PCA to the images to reduce data dimension, and then use SVM classifier to realize face recognition. Specifically, for a given face image, find who is in the image.

2 Tasks

Report the classification accuracy based on your training data by varying parameter K (i.e., the number of projections in PCA).

题目梗概

本次实验是一个关于人脸识别的实验，就是先要求你对于给的人脸的数据进行降维的处理，然后再通过SVM的方式进行区别。这边呢，我可以贴上一篇paper，是和我们这次实验有关的一篇paper，感觉还是写得不错的。

大致思想

一般来说，pca的步骤都大体相似，先去中心化，然后构造协方差矩阵，然后对这个协方差矩阵进行特征向量和特征值的分解。通过特征值的选取，比如最大的top100，这样就可以选取top100特征值对应的特征向量。然后我们有原数据比如200x10000的一个矩阵，10000个feature，200个sample，我们将这个10000x100的特征向量与原数据相乘200x10000*10000x100，即可将这个矩阵降维至100个feature。

实验

读取数据

读取pgm

clear;clc;
im1 = imread('./att_faces/orl_faces/s1/2.pgm');
imshow(im1);
m=im2double(im1);

这样我们就可以读取出112x92的一个格式为uint8的矩阵。

不过一般在matlab中,matrix的函数基本都是基于double类型的，所以我们可以对uint8类型的矩阵转换乘double，只需要使用函数im2double即可。

这边我实现的思路是，将训练集降维之后，本来是一个高维矩阵，我把它降到他自身的1/4。

然后得到一个训练集，现在就可以测试了。

如何区分

test读进来的肯定是一个原来维度的矩阵，还是先降维，降完维之后两种方式：欧式距离或者SVM进行区分

欧式距离

很简单，比如你有一个100个人的人脸训练数据，降完维度之后，你用test的降维后的矩阵和各个人的训练的矩阵求欧式距离，哪个最小，就是哪个人。

SVM

这个也简单，不过只能实现两个人之间的互判，就是说A的数据和B的数据分别打上1和-1的label，然后建立SVM，求max_{alpha}margin，然后根据alpha算出w和b,这样就能建立一个超平面，根据test的feature往里代入判断正负就ok。直接线性的SVM就ok，不用kernel升维度，已经实现了，但是代码中的一些细节部分还在整理，过几天附上来。

code

欧式距离法：

clear;clc;
PathRoot='./att_faces/orl_faces/s3/';
list=dir(fullfile(PathRoot));
fileNum=size(list,1)-3;
storematrix=zeros(56*46,fileNum);
sum0=zeros(56*46,1);
for k=3:fileNum+2
    path=list(k).name; % 这就是文件名，如果有子文件夹，则也包含在里面。
    path=[PathRoot,path];
    im1=imread(path);
    im1=imresize(im1,0.5);
    m=im2double(im1);
    B=reshape(m,56*46,1);
    sum0=B+sum0;
    storematrix(:,k-2)=B;
end
sum0=sum0./fileNum;
convex=zeros(56*46,56*46);
for i=1:fileNum
   storematrix(:,i)=storematrix(:,i)-sum0; 
end
convex=storematrix*storematrix';
[Vnew,Dnew]=eig(convex);
length0=length(Dnew);
judgemax=zeros(1,length0);
for i=1:length0
    judgemax(i)=Dnew(i,i);
end
[AS,pos]=sort(judgemax(:));
maxd=AS(3/4*length0:length0,1);
maxv=Vnew(:,pos(3/4*length0:length0,1));
V=storematrix'*maxv;
Vf=mean(V,1);
Vf=uint8(Vf);
% 欧式距离判别法
PathRoot1='./att_faces/orl_faces/s3/9.pgm';
imtest1=imread(PathRoot1);
imtest1=imresize(imtest1,0.5);
testm1=im2double(imtest1);
Btest1=reshape(testm1,56*46,1);
Btest1=Btest1-sum0;
Btest1Convex=Btest1*Btest1';
[Vnewt1,Dnewt1]=eig(Btest1Convex);
maxvt1=Vnewt1(:,pos(3/4*length0:length0,1));
vftest1=Btest1'*maxv;
count1=V;
[m1,n1]=size(V);
length1=m1;
for i=1:length1
    count1(i,:)=count1(i,:)-vftest1;
end
count1=count1.^2;
score1=sum(sum(count1));
PathRoot2='./att_faces/orl_faces/s5/3.pgm';
imtest2=imread(PathRoot2);
imtest2=imresize(imtest2,0.5);
testm2=im2double(imtest2);
Btest2=reshape(testm2,56*46,1);
Btest2=Btest2-sum0;
Btest2Convex=Btest2*Btest2';
[Vnewt2,Dnewt2]=eig(Btest2Convex);
maxvt2=Vnewt2(:,pos(3/4*length0:length0,1));
vftest2=Btest2'*maxv;
count2=V;
[m2,n2]=size(V);
length2=m2;
for i=1:length2
    count2(i,:)=count2(i,:)-vftest2;
end
count2=count2.^2;
score2=sum(sum(count2));

if score1<score2
    disp('判断出第一个测试者是训练集合中的人'); 
else
    disp('判断出第二个测试者是训练集合中的人');
end

SVM实现法，多重投票机制：

clear;clc;
lengthtrain=10;
ratio=0;
K=10;
sum0store=zeros(lengthtrain,lengthtrain,56*46);
omegastore=zeros(lengthtrain,lengthtrain,K);
bstore=zeros(lengthtrain,lengthtrain);
posstore=zeros(lengthtrain,lengthtrain,K);
maxvstore=zeros(lengthtrain,lengthtrain,56*46,K);
for train1=1:lengthtrain
train1str=num2str(train1);
pathr='./att_faces/orl_faces/s';
xiegang='/';
pathr=[pathr,train1str];
pathr=[pathr,xiegang];
PathRoot=pathr;
list=dir(fullfile(PathRoot));
fileNum=size(list,1)-7;% -2是读10个
for train_1=train1:lengthtrain
    if train_1==train1
       continue; 
    end
train_1str=num2str(train_1);
pathr_1='./att_faces/orl_faces/s';
xiegang='/';
pathr_1=[pathr_1,train_1str];
pathr_1=[pathr_1,xiegang];
PathRoot_2=pathr_1;
list_2=dir(fullfile(PathRoot_2));
fileNum_2=size(list_2,1)-7;% -2是读10个
storematrix=zeros(56*46,fileNum+fileNum_2);
sum0=zeros(56*46,1);

for k=3:fileNum+fileNum_2+2
    if k<=fileNum+2
        path=list(k).name; % 这就是文件名，如果有子文件夹，则也包含在里面。
        path=[PathRoot,path];
        im1=imread(path);
        im1=imresize(im1,0.5);
        m=im2double(im1);
        B=reshape(m,56*46,1);
        sum0=B+sum0;
        storematrix(:,k-2)=B;
    else
        path_2=list_2(k-fileNum).name;
        path_2=[PathRoot_2,path_2];
        im1=imread(path_2);
        im1=imresize(im1,0.5);
        m=im2double(im1);
        B=reshape(m,56*46,1);
        sum0=B+sum0;
        storematrix(:,k-2)=B;        
    end
end

sum0=sum0./(fileNum+fileNum_2);
convex=zeros(56*46,56*46);
for i=1:(fileNum+fileNum_2)
   storematrix(:,i)=storematrix(:,i)-sum0; 
end
sum0store(train1,train_1,:)=sum0;
convex=storematrix*storematrix';
[Vnew,Dnew]=eigs(convex,K);
length0=length(Dnew);
judgemax=zeros(1,length0);
for i=1:length0
    judgemax(i)=Dnew(i,i);
end
[AS,pos]=sort(judgemax(:));
maxd=AS(1:length0,1);
maxv=Vnew(:,pos(1:length0,1));
maxvstore(train1,train_1,:,:)=maxv;
V=storematrix'*maxv;
Vtrain1=V(1:fileNum,:);
% 训练集合1
Vtrain2=V(fileNum+1:fileNum+fileNum_2,:);
% 训练集合2

%先解出alpha,要想使用quadprog，得先将函数转化为对应形式
train=V;
m=fileNum+fileNum_2;
y=ones(m,1);
for i=fileNum+1:m 
    y(i,1)=-1;
end
f=zeros(1,m);
for i=1:m
f(1,i)=-1;
end
lengthx=length(train);
xlist=train;
H=zeros(m,m);
for i=1:m
for j=1:m
xi=zeros(1,lengthx);
xi(1,:)=xlist(i,:);
xit=xi';
xj=zeros(1,lengthx);
xj(1,:)=xlist(j,:);
H(i,j)=y(j,1)*y(i,1)*(xj*xit);
end
end
C=1000;
Aeq=y';
beq=0;
lb=zeros(m,1);
ub=C*ones(m,1);
alpha=quadprog(H,f,[],[],Aeq,beq,lb,ub);

omega=alpha.*y.*xlist;
omegaf=zeros(lengthx,1);
for i=1:lengthx
   omegaf(i,1)=sum(omega(:,i)); 
end
active=find(alpha>0);
sumact=length(active);
storey=zeros(sumact,1);
storex=zeros(sumact,lengthx);
for i=1:sumact
    storey(i,1)=y(active(i,1),1);
    storex(i,:)=train(active(i,1),:);
end
omegaf=omegaf';
storex=storex';
res=omegaf*storex;
res=storey-res';
fenzi=sum(res(:,1));
b=fenzi/sumact;
% SVM构建完毕
omegastore(train1,train_1,:)=omegaf;
bstore(train1,train_1,:)=b;
posstore(train1,train_1,:)=pos;
end
end

for testdata=1:lengthtrain
    mscore=zeros(1,lengthtrain);
    testdatastr=num2str(testdata);
    pathtest='./att_faces/orl_faces/s';
    whichone='/9.pgm';
    pathtest=[pathtest,testdatastr];
    pathtest=[pathtest,whichone];
    PathRoot1=pathtest;% test的数据
    imtest1=imread(PathRoot1);
    imtest1=imresize(imtest1,0.5);
    testm1=im2double(imtest1);
    Btest1=reshape(testm1,56*46,1);
    Btest1=Btest1-sum0;
    Btest1Convex=Btest1*Btest1';
    [Vnewt1,Dnewt1]=eigs(Btest1Convex,K);
    length0=length(Vnewt1);
    for t1=1:lengthtrain
        for t_1=t1:lengthtrain
            if(t1==t_1)
                continue;
            end
            Btest1=reshape(testm1,56*46,1);
            sum0=reshape(sum0store(t1,t_1,:),56*46,1);
            Btest1=Btest1-sum0;
            Btest1Convex=Btest1*Btest1';
            [Vnewt1,Dnewt1]=eigs(Btest1Convex,K);
            length0=K;
            pos=reshape(posstore(t1,t_1,:),1,K);
            maxvt1=Vnewt1(:,pos(1:length0));
            maxv=reshape(maxvstore(t1,t_1,:,:),56*46,K);
            vftest1=Btest1'*maxv;
            omegaf=reshape(omegastore(t1,t_1,:),1,K);
            stem=omegaf*vftest1';
            score=omegaf*vftest1'+bstore(t1,t_1);
            if score>0
               mscore(1,t1)=mscore(1,t1)+1; 
            else
               mscore(1,t_1)=mscore(1,t_1)+1;
            end
        end
    end
    [ansx,ansy]=find(mscore==max(max(mscore)));
    ansstr=['识别出是第',num2str(ansy)];
    ansstr=[ansstr,'个人'];
    disp(ansstr);
    if ansy==testdata
       ratio=ratio+1; 
    end

end
ratio=ratio/lengthtrain;
ratiostr=['准确率为',num2str(ratio)];
disp(ratiostr);

转载于:https://www.cnblogs.com/harrysong666/p/10170744.html