ipca_jc_1.m

一种增量的人脸识别算法——增量PCA学习算法（matlab实现）

% IPCA_JC_1.m   19th, May, 2008  made by J.Chen
% reference: 一种增量PCA算法及其在人脸识别中的应用. 夏鹏，张浩然，徐展敏

%-------------------------------PCA & Test---------------------------------

X1 = [];                                                                   %所有训练图片

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:200
    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,namda,V] = svd(X1); % put in 1/sqrt(n-1) later
dnamda = diag(namda);

%以下选择90%的能量--------------
dsum = 0;
for i=1:200  
    dsum = dsum + dnamda(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 + dnamda(i)^2;
  end
end

%选择k个主分量------------------
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=6: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:200
            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; 
   end; 
end; 

accuracy = accu/200                                                        % 输出识别率

%-------------------------------PCA & Test---------------------------------

%-------------------------------IPCA---------------------------------------

X2 = [];                                                                   %新增的训练图片r个

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']';                                                    %X2是一个N*r矩阵，X2中每一列数据代表一张图片，其中r=40,size(X2)=[10304 40]
    end
end

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

rnamda = zeros(k,k);
for i=1:k
    rnamda(i,i) = dnamda(i).^(-0.5);
end
H1 = U*rnamda;                                                             %size(H1)=[n k]

%----------------------------------
%S2 = X2*X2';                                                               %size(S2)=[n n]
%S2bar = H1'*S2*H1;                                                         %size(S2bar)=[k k]

A = H1'*X2;
B = A*X2';
S2bar = B*H1;
%----------------------------------

[P2 namda2] = eig(S2bar);                                                  %size(P2)=[k k]
dnamda2 = diag(namda2);
P = H1*P2;                                                                 %size(P)=[n k]  找到了S1+S2的特征向量P，即更新后的特征空间。

X = [];
Y = [];
X = [X1';X2']';
Y = P'*X;

%测试过程---------------------------

accu2 = 0;

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

accuracy2 = accu2/160                                                      % 输出识别率

%-------------------------------IPCA---------------------------------------