pca_jc.m

经典人脸识别算法——PCA（Matlab实现）

% PCA_JC.m   18th,May,2008  made by J.Chen

X1 = [];                                                                     %所有训练图片
m = 200;
for i = 1:40
    for j = 1:5
        a = imread(strcat('E:\FaceDataBase\ORL\s',num2str(i),'\',num2str(j),'.pgm'));
                                                                           %a是112*92的矩阵
        b = (a(1:112*92))';                                                   %b是列向量N*1，N=10304，提取顺序是先列后行，
                                                                           %即从上到下，从左到右
        b = double(b);
        X1 = [X1';b']';                                                       %X1是一个N*M矩阵，X1中每一列数据代表一张图片，其中M=200,size(X1)=[10304 200]
    end
end

%中心化数据
Xm = mean(X1');                                                                  %X1转置后每行代表一张图片，按列求平均,返回一个行向量1*10304
for j = 1:m
    Xcter(:,j) = X1(:,j)-Xm';
end
X1=Xcter;                              %此时，X1为中心化的样本数据

% The principal component coefficients are the eigenvectors of
% S = X1'*X1./(n-1), but computed using SVD.
[U0,sigma,V] = svd(X1); % put in 1/sqrt(n-1) later
dsigma = diag(sigma);


%dsigma=flipud(dsigma);
%U0=fliplr(U0);


%以下选择90%的能量
dsum = 0;
for i=1:m  
    dsum = dsum + dsigma(i)^2;
end
dsum_extract=0;
k=0;
while(dsum_extract/dsum<0.9)
  dsum_extract = 0;
  k=k+1;
  for i=1:k  
      dsum_extract = dsum_extract + dsigma(i)^2;
  end
end                                    %选择k个主分量

%%%%%%k = 200;
for i=1:k
    U(:,i) = U0(:,i);
end

Y = U'*X1;                             %size(Y)=[k M],是图片向量映射到k维特征空间的坐标

accu = 0;                                            %下面的人脸识别过程中就是利用这些组合系数来进行识别
 
%测试过程

for i=1:40
    for j=9:10                             %读入40 x 5 副测试图像
        a=imread(strcat('E:\FaceDataBase\ORL\s',num2str(i),'\',num2str(j),'.pgm'));
        b=(a(1:10304))';
        b=double(b);
        bm = b - Xm';                       %将测试图片中心化
        t=U'*bm;                            %计算测试图片在新特征空间中的坐标，是k*1阶矩阵
        for p=1:m
            mdist(p)=norm(t'-Y(:,p)');
        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; 
      if class1==i
          accu = accu+1;
      end
   end; 
end; 

accuracy = accu/80                  % 输出识别率

%下面进行增量IPCA

%for i = 1:k
%    for j=1:k
%        ssigma(i,j) = sigma(i,j).^2;
%    end        
%end

%for i=1:k
%    VV(:,i) = V(:,i);
%end
%W = X1*VV;
%H1 = W.*((200)^(1/2))*(ssigma.^(-1));
%X2 = [];                                                                     %所有训练图片

%for i = 1:40   %增量
%    for j = 6:6
%        a = imread(strcat('E:\FaceDataBase\ORL\s',num2str(i),'\',num2str(j),'.pgm'));
                                                                           %a是112*92的矩阵
%        b = (a(1:112*92))';                                                   %b是列向量N*1，N=10304，提取顺序是先列后行，
                                                                           %即从上到下，从左到右
%        b = double(b);
%        X2 = [X2';b']';                                                       %X1是一个N*M矩阵，X1中每一列数据代表一张图片，其中M=200,size(X1)=[10304 200]
%    end
%end

%中心化数据
%Xm = mean(X2');                                                                  %X1转置后每行代表一张图片，按列求平均,返回一个行向量1*10304
%for j = 1:40   %增量
%    Xcter(:,j) = X2(:,j)-Xm';
%end
%X2 = Xcter;                              %此时，X2为中心化的样本数据