featureextract.m

这是人脸识别代码

function [W,ProjectCoef,ProjectCoefClass,ExtractTime]=FeatureExtract(X,XClass,Judge,sDatabase,sOrder)
%
% [W,ProjectCoef,ProjectCoefClass,ExtractTime]=FeatureExtract(X,XClass,Judge)
%
% 输入样本按列堆积的矩阵X，以及每一列的所属的类别向量XClass，根据判据Judge提取特征。
% 当Judge='PCA'时（默认），求得并选取p个总体散度矩阵St的特征向量组成Eigen脸W，使得保留90%的能量。
% 当Judge='LDA'时，求得并选取L=(K-1)个(Sw^-1)*Sb的特征向量组成Fisher脸W， K是类别数。
% Sb是样本X的类间散度矩阵，Sw是样本X的类内散度矩阵。XClass(i)是样本X(:,i)的类别。
% 并返回每类样本在W上的投影系数ProjectCoef和所花费的时间ExtractTime。
% ProjectCoef(:,i)所属的类别为ProjectCoefClass(i)。
% 返回该投影系数是因为一般来说会用每类样本各自取平均在W上投影作为一类样本的特征。
% 
% Author : kk.h
% Date : 2004.5.2
% Email : [e-mail]
% SUN YAT-SEN UNIVERSITY 
%

% -----------------------------------------------------------------------------------------------------------

% 开始计算时间
BTime = clock;

% -----------------------------------------------------------------------------------------------------------
% 得到相应的.mat文件
ssDatabase = sDatabase;
ssDatabase = strrep(ssDatabase,':','_');
ssDatabase = strrep(ssDatabase,'\','_');

sMatFile = fullfile('mat',[ssDatabase '_' sOrder '_train_' Judge '.mat']);

% 如果相应的.mat文件存在，直接load该文件，得到变量返回。不必重新读取
dirMatFile = dir(sMatFile);
if size(dirMatFile,1)~=0
    load(sMatFile);
    ExtractTime = etime(clock,BTime);
    return;
end;

% -----------------------------------------------------------------------------------------------------------
% 样本数
XCount = size(X,2); 

% 计算每一类的均值。
[XMean,XMeanClass]= MeanClass(X,XClass);

% % 初始化
% ClassCount = 1; % 类别数
% XMeanClass(1) = XClass(1); % 第1类的类别
% XMean(:,1) = X(:,1); % 第1类的样本总和
% XMeanCount(1) = 1; % 第1类的样本数量
% 
% for i=2:XCount
%     if length(strfind(XMeanClass,XClass(i))) == 0 % 如果第i个样本的类别不在已有的类别中
%         ClassCount = ClassCount + 1;
%         XMeanClass(ClassCount) = XClass(i); % 第ClassCount类的类别
%         XMean(:,ClassCount) = X(:,i); % 第ClassCount类的样本总和
%         XMeanCount(ClassCount) = 1; % 第ClassCount类的样本数量
%     else
%         XMean(:,ClassCount) = XMean(:,ClassCount) + X(:,i);  % 第ClassCount类的样本总和
%         XMeanCount(ClassCount) = XMeanCount(ClassCount) + 1;  % 第ClassCount类的样本数量
%     end
% end;
% 
% % 每一类分别除去各自的数量得到每一类的均值
% for i=1:ClassCount
%     XMean(:,i) = XMean(:,i) / XMeanCount(i);
% end;

% -----------------------------------------------------------------------------------------------------------
% 总体均值
m = mean(X,2);

% 总体散度矩阵： nxn 维
Xt = X - kron(m,ones(1,XCount));
% St = Xt * Xt';  % St = Sw + Sb;

% 每一样本减去各自所属类的均值。
for i=1:XCount
    Xw(:,i) = X(:,i) - XMean(:,strfind(XMeanClass,XClass(i)));
end;

% 类内散度矩阵： nxN * Nxn = nxn 维。n是每个样本的大小，N是样本总数  rank(Sw) <= N-K
% Sw = Xw * Xw';

% 每一类的均值分别减去总体均值
for i=1:ClassCount
    Xb(:,i) = XMean(:,i) - m;
end;

% 类间散度矩阵： nxK * Kxn = nxn 维。n是每个样本的大小，K是类别数。 rank(Sb) <= K-1
% Sb = Xb * Xb' ;
% -----------------------------------------------------------------------------------------------------------
% 求解总体散度矩阵St或者类间散度矩阵Sb的特征向量，用来降维压缩数据
% 但 Sc=(1/N)Xc*Xc' 的维数太大，根据SVD定理的推论，可以先求较低维的 Rc=Xc'*Xc 的特征向量组V。
% ？？(1/N)Xc*Xc'和Xc*Xc'的正交归一的特征向量EigenFace和U是一样的吗？
% Xc*Xc'是对称半正定的矩阵，不会有负的特征根

Xc = Xt;  % Xc：nxN 用总体散度矩阵St作产生矩阵
% Xc = Xb;  % Xc：nxK 用类间散度矩阵Sb作产生矩阵

% Rc的维数： (nxN)' * nxN = NxN （n是样本的维数，N是样本总数XCount）
% Rc的维数： (nxK)' * nxK = KxK （n是样本的维数，K是样本类别数）
Rc = Xc'*Xc; 

% 求解Rc=Xc'*Xc的特征根Vc
[Vc,Dc] = eig(Rc); % Rc*Vc = Vc*Dc

% 直接做SVD内存不足
% [Uc,Dc,Vc] = svd(Xc);

% -----------------------------------------------------------------------------------------------------------
% 按St的特征值从大到小排序

% 对角向量
DcCol = diag(Dc);

% 从小到大排序
[DcSort,DcIndex] = sort(DcCol);

% Vc的列数
VcCols = size(Vc,2);

% 反序
for (i=1:VcCols)
    VcSort(:,i) = Vc(:,DcIndex(VcCols-i+1));
    DcSort(i) = DcCol(DcIndex(VcCols-i+1));
end

% -----------------------------------------------------------------------------------------------------------
% 降维

% 如果是LDA，则要用PCA先降维，即用总体散度矩阵的最大的前p个特征向量组成Eigen脸W，使得|W'StW|最大。
% W的维数： n x p ，其中： p = 样本个数N - 类别数K
if Judge=='LDA'
    p = XCount-ClassCount;
 %   p = ClassCount; % 用类间散度矩阵Sb来降维，rank(Sb) <= K-1。
else % 如果是PCA，则默认保留90%的能量
    DcSum = sum(DcSort);
    DcSum_extract = 0;
    p = 0;

    while( DcSum_extract/DcSum < 0.9)
        p = p + 1;
        DcSum_extract = sum(DcSort(1:p));
    end
end;

% -----------------------------------------------------------------------------------------------------------
% Wpca是由前p个最大的非0特征值对应的特征向量组成的，根据SVD定理从Vt计算得到。维数为 nxp
i=1;
while (i<=p && DcSort(i)>0)
    Wpca(:,i) = DcSort(i)^(-1/2) * Xc * VcSort(:,i); 
    i = i + 1;
end

% -----------------------------------------------------------------------------------------------------------
% -----------------------------------------------------------------------------------------------------------
% 直接做LDA会出现小样本问题，即当 n<N时，即样本数小于样本维数时，Sw是奇异的，无法求解(Sw^-1)*Sb的特征向量。
% 所以先用PCA降维。 Sw，Sb 都由 nxn 降为 pxp
if Judge=='LDA'
    Sw = Wpca' *  Xw * Xw' * Wpca; % (nxp)' * nxN * Nxn * nxp = pxp
    Sb = Wpca' *  Xb * Xb' * Wpca; % (nxp)' * nxK * Kxn * nxp = pxp

    % 按Fisher法则求解pxp维矩阵的特征向量
    [Vw,Dw] = eig((Sw^-1)*Sb); 

    % 按特征值从大到小排序

    % 对角向量
    DwCol = diag(Dw);

    % 从小到大排序
    [DwSort,DwIndex] = sort(DwCol);

    % 反序
    for (i=1:p)
        VwSort(:,i) = Vw(:,DwIndex(p-i+1));
        DwSort(i) = DwCol(DwIndex(p-i+1));
    end

    % 降维，保持99%的能量依然只是从160个中取到10个。这是因为这时候特征值的大小已经和能量没什么关系。
    % 但是由于 rank((Sw^-1)*Sb) <= rank(Sb) <= K-1，所以非零的特征根的个数最多为K-1。所以可以到 K-1 维。
    % 其它特征根都是0，对应的特征向量对类间已经没有什么区分作用。
    L = ClassCount - 1;
    
    Wlda = VwSort(:,1:L);  % 维数： pxL （即行是PCA降维后的大小，列是LDA降维之后的特征向量个数L）
    
    W = Wpca * Wlda;  % 维数： nxp * pxL = n x L （即列是样本维数，行是LDA降维后的特征向量个数L(<=K-1)）
    
else
    W = Wpca;
end;
% -----------------------------------------------------------------------------------------------------------
% 结果数据

% % 所有训练样本都在特征脸上投影，每一列的系数代表每一个样本的特征
ProjectCoef = W'*X;  % 维数： (nxL)' * nxN = LxN
ProjectCoefClass = XClass;

% 用每一类的均值在特征脸上投影，投影每一列的系数代表每一类样本的特征
% ProjectCoef = W'*XMean;  % 维数： (nxL)' * nxK = LxK
% ProjectCoefClass = XMeanClass;

% 变换矩阵W=Wpca*Wlda把总体离差阵也变成对角阵了。因为Wpca首先把St变成对角的了。所以这样做出来的LDA也是统计不相关的。

% 所用总时间
ExtractTime = etime(clock,BTime);


% -----------------------------------------------------------------------------------------------------------
% 保存结果。

% 如果mat的文件夹不存在
if size(dir('mat'),1)==0
  mkdir('mat');
end;

save(sMatFile,'W','ProjectCoef','ProjectCoefClass');