orl_lda.m

用MATLAB编写的人脸检测算法

clear all;
close all;
ratio = 2;
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_2;

disp('orl fisher face now begin');
right_num = [];
for count = 5 %2:8 
    count
	train_num = count;
    train_total = train_num * class_num;
	test_num = sample_num - train_num;
    test_total = test_num * class_num;
    
    %% -------- 1.calculate Sb Sw -------- %%
    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));
            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));
            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        
	
    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;
    
    %% -------- 2.calculate transform of St -------- %%
    ss = (St' * St) / train_total;   %% [train_total, train_total]
%    rr = rank(ss);
    rr = rank(ss)-class_num;
    
    [V_St, D_St] = eig(ss);
    dd = abs(eig(ss));
    [dd_value, dd_site] = sort(dd);
    temp2 = rr;
    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;
    
    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;

    %% -------- 3.calculate transform of unified fisher vectors -------- %%
    ss = inv(Sww)*Sbb;

    [V_ss, D_ss] = eig(ss);
    dd = abs(eig(ss));
    [dd_value, dd_site] = sort(dd);
    dd_value';
    temp1 = rr - (class_num-1) + 1;
    yy = [];
    for i = rr:-1:temp1
        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)';

    for ii = 1:10
        ww_temp = ww(ii,:);
        ww_temp = ww_temp - min(ww_temp);
        ww_temp = ww_temp / max(ww_temp);
        ww_temp = floor(ww_temp * 255);
        orl_temp = zeros(row, line);
        for p = 1:row
            for q = 1:line
                  site = (p-1)*line + q;  
                  orl_temp(p,q) = ww_temp(site);
            end
        end
        orl_temp = uint8(orl_temp);
        ss = sprintf('fisherorl_%d.bmp',ii);
        imwrite(orl_temp, ss,'bmp');
    end    
    disp('orl_fisherface is over');
    
    %% -------- 4.calculate transform of test and train -------- %%
    orl_train_new = ww * orl_train;
    orl_test_new = ww * orl_test;

	error = [];
	for i = 1:class_num
        for j = 1:test_num
            %i,j
            kk = (i-1)*test_num + j;
            res_y = orl_test_new(:,kk);
        
            dis = [];
            for ii = 1:class_num
                for jj = 1:train_num
                    res_x = orl_train_new(:,(ii-1)*train_num+jj);
                    res = (res_y - res_x).^2;
                    temp = sqrt(sum(res));
                    dis = [dis; temp];
                end
            end
            [value, site] = min(dis);
            site_new = floor((site-1)/train_num) + 1;
            
            if site_new ~= i
                error = [error; i*test_num+j];
            end
        end
	end
    right_num = [right_num; (1.0-length(error)/(class_num*test_num))*100 ];
end

right_num'