orl_rbf_0822_1534.m

MATLAB做的不错的人脸识别的程序

%% derived from "Face Recognition With Radial Basis Function (RBF) Neural Networks" 
%% which is written by M.J.Er,S.Q.Wu, J.W.Lu and H.L.Toh in 2002. 

clear all;
close all;
ratio = 4;
class_num = 40;
sample_num = 10;
row_0 = 112;
line_0 = 92;
row = ceil(row_0 / ratio);
line = ceil(line_0 / ratio);
row_line = min(row, line);
rowline = row * line;
load orl_4;

right_num_test = [];
right_num_train = [];

disp('orl fisherface now begin');
right_num = [];
for count = 5
    tic;
    
    count
	train_num = count;
    train_total = train_num * class_num;
	test_num = sample_num - train_num;
    test_total = test_num * class_num;
    
    
    orl_train = zeros(rowline, train_total);
    orl_m_all = zeros(rowline,1);
    orl_m = zeros(rowline,class_num);
	for i = 1:class_num
        for j = 1:train_num
            kk = (i-1)*train_num + j;
            orl_temp = double(orl(:,:,j,i));
            %orl_temp = double(orl(:,:,j+test_num,i));
            for p = 1:row
                for q = 1:line
                      site = (p-1)*line + q;  
                      orl_train(site, kk) = orl_temp(p,q);
                end
            end
            orl_m_all = orl_m_all + orl_train(:,kk);
            orl_m(:,i) = orl_m(:,i) + orl_train(:,kk);
        end
        orl_m(:,i) = orl_m(:,i)/train_num;
    end
    orl_m_all = orl_m_all / train_total;
    
    orl_test = zeros(rowline, test_total);
  	for i = 1:class_num
        for j = 1:test_num
            kk = (i-1)*test_num + j;
            orl_temp = double(orl(:,:,j+train_num,i));
            %orl_temp = double(orl(:,:,j,i));
            for p = 1:row
                for q = 1:line
                    site = (p-1)*line + q;  
                    orl_test(site, kk) = orl_temp(p,q);
                end
            end
        end
    end        
	%% -------- 1.Perform fisherface method to reduce the demension of the feature space -------- %%
        %%------- PCA -------%%
    Sb = zeros(rowline,class_num);
    for i = 1:class_num
        Sb(:,i) = orl_m(:,i) - orl_m_all;
    end
%    Sb0 = (Sb * Sb') / class_num;  %% [rowline, rowline]      

    Sw = zeros(rowline,train_total);
	for i = 1:class_num
        for j = 1:train_num
            kk = (i-1)*train_num + j;
            Sw(:,kk) = orl_train(:,kk) - orl_m(:,i);
        end    
    end

    St = zeros(rowline,train_total);
	for i = 1:class_num
        for j = 1:train_num
            kk = (i-1)*train_num + j;
            St(:,kk) = orl_train(:,kk) - orl_m_all;
        end    
    end
%    St0 = (St * St') / train_total;  %% [rowline, rowline]
    
    clear orl_m_all orl_m;
    
    ss1 = (St' * St) / train_total;   %% [train_total, train_total]
%    rr = rank(ss);
    for rr1 = 164%rank(ss1)-1*class_num%:rank(ss1);%rank(ss1);%46;%
    
    [V_St, D_St] = eig(ss1);
    dd = abs(eig(ss1));
    [dd_value, dd_site] = sort(dd);
    temp2 = rr1;
    temp1 = train_total - temp2 + 1;
    yy = [];
    Dt = zeros(temp2);
    jj = 0;
    for i = train_total:-1:temp1
        yy = [yy; V_St(:,dd_site(i))'];
        jj = jj + 1;
        Dt(jj,jj) = dd_value(i)^(-0.5);   
    end
    yy = yy';
    zz = St * yy * Dt;
    disp('PCA is over');
    
    Stt = zz' * St * St' / train_total * zz / train_total;
    Sww = zz' * Sw * Sw' / train_total * zz / train_total;    
    Sbb = zz' * Sb * Sb' / class_num * zz / train_total;

        %% -------- LDA -------- %%
    ss2 = inv(Stt)*Sbb;
%    ss = inv(Sww)*Sbb;

    [V_ss, D_ss] = eig(ss2);
    dd = abs(eig(ss2));
    [dd_value, dd_site] = sort(dd);
    for temp3 = class_num - 1; %rank(ss)
    temp4 = rr1 - temp3 + 1;
    
    yy = [];
    for i = rr1:-1:temp4
        v_temp = V_ss(:,dd_site(i))';
        v_temp = v_temp / norm(v_temp);
        yy = [yy; v_temp];
    end
    yy = yy';
    v0 = yy;
    ww = (zz * v0)';
    
    dd_bak = dd_value(temp1:rr1)';
    disp('LDA is over');
    toc;
    
    y_train = ww * orl_train;
    y_test = ww * orl_test;
    %% -------- 2. Perform RBF networks -------- %%
    tic;
    u = class_num; % the number of RBF units
    clus = zeros(1,train_total);
    for i = 1:train_total
        kk = floor((i-1)/train_num) + 1;
        clus(i) = kk;
    end
    flg = zeros(1,u);
    FlgMin = min(flg);      
    Epoch1 = 1;
    while (FlgMin == 0) & (Epoch1<=100)
          Epoch1 
          flg = zeros(1,u);
          y_m = zeros(class_num-1,u);
          num = ones(1,u);
          for i = 2:length(clus)
              if clus(i) == clus(i-1)
                  num(clus(i-1)) = num(clus(i-1))+1;
              end
          end
          kk = 0;
          for i = 1:u
              for j = 1:num(i)
                  kk = kk +1;
                  y_m(:,i) = y_m(:,i) + y_train(:,kk);
              end
              y_m(:,i) = y_m(:,i)/num(i);
          end
          Dis = zeros(1,u);% to find the distance between y_m(:,i) and the furthest ponit from y_m(:,i) in class i  
          kk = 0;
          for i = 1:u
              for j = 1:num(i)
                  kk = kk+1;
                  d_temp = sqrt((y_train(:,kk)-y_m(:,i))'*(y_train(:,kk)-y_m(:,i)));
                  if d_temp >= Dis(i)
                      Dis(i) = d_temp;
                  end
              end
           end
           DisC = zeros(u,u);% distance between y_m(:,i) and y_m(:,j)
           DisI = zeros(1,u);% index of the nearest point from the mean of class i
           DisMin = zeros(1,u);
           uu = u;
           for i = 1:u
               DisMin(i) = sqrt((y_m(:,i)-y_m(:,u-i+1))'*(y_m(:,i)-y_m(:,u-i+1)));
               DisI(i) = u-i+1;
               for j = 1:u
                   DisC(i,j) = sqrt((y_m(:,i)-y_m(:,j))'*(y_m(:,i)-y_m(:,j)));
                   if DisC(i,j) <= DisMin(i) && DisC(i,j)~=0
                      DisMin(i) = DisC(i,j);
                      DisI(i) = j;
                  end
               end
               l = DisI(i);
               if (Dis(i)+Dis(l))<= DisMin(i)
                   flg(i) = 1;
               end 
               if ((Dis(i)+Dis(l))>DisMin(i)) && ((Dis(i)-Dis(l))>0)
                   P = randperm(4); k = P(1); %% choose a number between 1 and 4 randomly
                   uu = uu+1;
                   point = sum(num(1:i-1))+k+1;
                   for jj = train_total:-1:point
                       clus(jj) = clus(jj)+1;
                   end
               end
           end
           u = uu;
           FlgMin = min(flg);
           Epoch1 = Epoch1+1;
    end
    Center = y_m; [M N] = size(Center);
    beta = 0.5 ; 
        eta = 0.7  ;
    DelW = zeros(1,u); DelB = zeros(1,u); Del = zeros(1,u);
    for i = 1:u
        DelW(i) = Dis(i)/sqrt(abs(log(beta)));
        DelB(i) = eta*DisMin(i);
        if DelB(i)> DelW(i)
            Del(i) = DelB(i);
        else
            Del(i) = DelW(i);
        end
    end
    W = zeros(class_num,u);
    T = zeros(class_num,train_total);
    for i = 1:class_num
        for j = 1:train_num
            kk = (i-1)*train_num+j;
            T(i,kk) = i;
        end
    end
    E = ones(1,class_num);
    Epoch2 = 1;
    while (max(E) >= 10.^-2) && (Epoch2<=40)
        Epoch2 
        E = zeros(1,class_num);
        Rbf_train = zeros(u,train_total);
        for i = 1:train_total
            for j = 1:u
                Rbf_train(j,i) = exp(-((y_train(:,i) - Center(:,j))'*(y_train(:,i) - Center(:,j)))/(2*Del(j).^2));
            end 
        end
        W = T*pinv(Rbf_train);
        Y = W * Rbf_train;
        LernRate =2*10.^1;
        for i = 1:class_num
            for j = 1:train_num
                ll = (i-1)*train_num + j;
                E(i) = E(i) + (T(:,ll)-Y(:,ll))'*(T(:,ll)-Y(:,ll))/2;
                %W(:,i) = W(:,i) + (T(:,ll)-Y(:,ll))'*Rbf_train(:,ll);
                kk = 0;
                ErrRatCen = zeros(M,N);
                ErrRatDel = zeros(1,u);
                kk = 0;
                for ii = 1:u
                    for jj = 1:num(ii)
                        kk = kk+1;
                        ErrRatCen(:,ii) = (2*LernRate*(T(:,kk)-Y(:,kk))'*W(:,ii)*Rbf_train(ii,kk)*(y_train(:,kk)-Center(:,ii))'/(Del(ii).^2))';
                        ErrRatDel(ii) = (2*LernRate*(T(:,kk)-Y(:,kk))'*W(:,ii)*Rbf_train(ii,kk)*((y_train(:,kk)-Center(:,ii))'*(y_train(:,kk)-Center(:,ii)))/(Del(ii).^3));
                    end
                end
                Center = Center + ErrRatCen;
                Del = Del + ErrRatDel; 
             end
         end
         Epoch2 = Epoch2+1;
     end
     Rbf_test = zeros(u,test_total);
     for i = 1:test_total
         for j = 1:u
             Rbf_test(j,i) = exp(-((y_test(:,i) - Center(:,j))'*(y_test(:,i) - Center(:,j)))/(2*Del(j).^2));
         end 
     end
    toc;
    disp('Radial Basis Function is over');
    
    %% -------- 4.calculate transform of test and train -------- %%
    orl_train_new = W * Rbf_train;
    orl_test_new = W * Rbf_test;

	error_test = [];
	[Y_test I_test] = max(orl_test_new,[],1);
    for i = 1:test_total
        j = floor((i-1)/test_num)+1;
        if I_test(i) ~= j
            error_test = [error_test; (i-1)*test_num+j];
        end
    end
    right_num_test = [right_num_test; (1.0-length(error_test)/(class_num*test_num))*100 ];
    right_num_test'
    error_test';
    error_train = [];
	[Y_train I_train] = max(orl_train_new,[],1);
    for i = 1:train_total
        j = floor((i-1)/train_num)+1;
        if I_train(i) ~= j
            error_train = [error_train; (i-1)*train_num+j];
        end
    end
    right_num_train = [right_num_train; (1.0-length(error_train)/(class_num*train_num))*100 ];
    right_num_train'
end
end
end