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 : kkcocoon@163.com
% 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');