OpenCV人臉識別Eigen算法源碼分析

1 理論基礎

學習Eigen人臉識別算法需要了解一下它用到的幾個理論基礎，現總結如下：

1.1 協方差矩陣

 首先需要了解一下公式：

共公式可以看出：均值描述的是樣本集合的平均值，而標准差描述的則是樣本集合的各個樣本點到均值的距離之平均。以一個國家國民收入為例，均值反映了平均收入，而均方差/方差則反映了貧富差距，如果兩個國家國民收入均值相等，則標准差越大說明國家的國民收入越不均衡，貧富差距較大。以上公式都是用來描述一維數據量的，把方差公式推廣到二維，則可得到協方差公式：

協方差表明了兩個隨機變量之間的相關性，值為正說明兩者是正相關的，值為負說明兩者是負相關的，值為零說明兩者不相關，舉一個簡單的小例子，假設一個人用4個維度身高、體重、距離屋頂的高度、每天畫畫的時間來表示：身高取樣X=[1 2 3 4 5 6 7 8 9]，體重取樣Y=[11 12 13 14 15 16 17 18 19]，距離屋頂的高度取樣Z=[9 8 7 6 5 4 3 2 1]，每天畫畫時間L=[1 1 1 1 1 1 1 1 1]，則有cov(X,Y)=7.5，cov(X,Z)=-7.5，cov(X,L)=0，結果很明顯X和Y協方差為正數兩者正相關，X和Z協方差為負數兩者負相關，X和L協方差為0，說明它們不相關。以上例子每一個隨機變量都可以表示一個維度，我們計算了部分維度之間的協方差，計算所有維度之間的協方差並組織成矩陣的形式，就有了協方差矩陣的概念：C_nxn=[c_i,j]=[cov(Dim_i,Dim_j)]   i,j=1,2,…,n，Dim_i表示第i個維度向量。以Matlab協方差矩陣為例，將X，Y，Z，L分別作為1，2，3，4個維度，則有c_1,1=7.5，c_1,2=7.5，c_1,3=-7.5，c_1,4=7.5……，所以協方差矩陣為：

 

在Matlab中可以把矩陣的每行看做是4個隨機變量的一組取樣樣本，每列看做是一個維度，則可以直接用con函數求得4個維度的協方差矩陣：

1.2 Jacobi迭代法求對稱矩陣特征向量及特征值

 雅可比迭代法的基本思想是：通過一組平面旋轉變換（相似正交變換）化對稱矩陣A為對角矩陣，進而求出A的特征值與特征向量。由線性代數理論可知：若矩陣A是實對稱矩陣，則一定存在正交矩陣U，使得U^T*A*U=D，其中D對角矩陣，其主對角線元素λ_i是A的特征值，正交矩陣U的第i列是A對應特征值λ_i的特征向量。於是求對稱矩陣A的特征值問題轉化為尋找正交矩陣U，使得U^T*A*U為對角矩陣，這個問題的困難在於如何構造U，為此我們先看一下平面上的旋轉變換：

則有：

其中：

上述推導其實說明了一種構造正交矩陣P，並使得P^T*A*P為對角矩陣的方法，可以將這種方法推廣到nxn對角矩陣，首先引入n階旋轉矩陣（Givens矩陣）的概念：

平面旋轉矩陣有如下性質：

（1）U_pq為正交矩陣，即U_pq^T*U_pq=E

^（2）U^TAU=B仍為對稱矩陣，且B與A有相同的特征值

Jacobi迭代法，在每一次迭代時都是進行一次（2）中的轉換，這里p、q分別是前一次的迭代矩陣A的非主對角線上絕對值最大元素的行列號，變換后元素值可以由以下公式求出：

由公式可以看出轉換后矩陣相比原矩陣只是在p，q行和列的元素發生了改變，旋轉角的計算過程和2維時一樣，其意義是使得apq和aqp值為零，這樣每次迭代都使得非對角線上絕對值最大的元素變為零，所以整個迭代的過程就是使對角線外元素逐步逼近於零，這是對角線上的元素即為原對稱矩陣的特征值λ_i。在進行Jacobi迭代時，假如i次迭代時旋轉矩陣為Ui，每次迭代對單位矩陣I依次左乘Ui，最終迭代結束后可得矩陣D=U_k…U₂U₁I，這里k為迭代次數，則可以證明D的列向量即為特征值λ_i對應的特征向量，證明如下：

上述推導過程中d_i為矩陣D的i列表示的列向量，由最后的等式及特征值定義，可以得知λ_i是A的特征值，d_i為對應的特征向量。

2 OpenCV源碼解析

2.1 關鍵函數

（1）void reduce(InputArray src, OutputArray dst, int dim, int rtype, int dtype=-1)

其英文注釋：transforms 2D matrix to 1D row or column vector by taking sum, minimum, maximum or mean value over all the rows.

其英文注釋不太准確，函數的作用其實是：將2維矩陣轉換為1維行向量或列向量，如轉換為行向量，則每列處的值為原矩陣對應列所有值的和，最小值，最大值，平均值；如轉換為列向量，則每行處的值為原矩陣對應行所有值的和。該函數參數意義如下：

src: 原矩陣

dst: 目的向量

dim: 指明處理后向量是行向量還是列向量，0原矩陣被處理成行向量，否則原矩陣被處理成列向量

op: 取值為CV_REDUCE_SUM，CV_REDUCE_MAX，CV_REDUCE_MIN，CV_REDUCE_AVG之一

dtype: 目的向量類型

（2）void gemm(InputArray src1, InputArray src2, double alpha, InputArray src3, double gamma, OutputArray dst, int flags=0)

其英文注釋：implements generalized matrix product algorithm GEMM from BLAS.

函數的作用：實現廣義矩陣乘法，只對最后一個參數進行說明

flags: 取值為GEMM_1_T，GEMM_2_T，GEMM_3_T之1或者它們的組合，例如取值為GEMM_1_T則進行乘法之前對src1進行轉置，所有函數作用可由以下公式來說明：

dst=alpha*op(src1)*op(src2)+gamma*op(src3)，其中op(X)是X還是X^T由flags確定。

（3）void mulTransposed( InputArray src, OutputArray dst, bool aTa, InputArray delta=noArray(), double scale=1, int dtype=-1 )

其英文注釋：multiplies matrix by its transposition from the left or from the right.

函數的作用：矩陣左乘或右乘其轉置矩陣，參數意思如下：

src: 原矩陣

dst: 目的矩陣

ata: 乘法順序，true A^T*A false A*A^T

delta：在進行乘法前src先減去該數組

scale：乘法之后對結果進行scale倍縮放

dtype：目的矩陣類型

當ata為真時可用公式 dst=(src-delta)^T*(src-delta)*scale 來說明函數的作用，該函數內部調用了函數（2）

（4）void calcCovarMatrix( InputArray samples, OutputArray covar, OutputArray mean, int flags, int ctype=CV_64F)

其英文注釋：computes covariation matrix of a set of samples

函數作用：計算矩陣行向量或列向量的協方差矩陣，該函數中會調用函數（3）來實現相應功能

（5）bool eigen(InputArray src, OutputArray eigenvalues, OutputArray eigenvectors, int lowindex=-1, int highindex=-1)

其英文解釋：finds eigenvalues and eigenvectors of a symmetric matrix

函數作用：求對稱矩陣的特征值和特征向量，在該函數中會利用Jacobi方法來求對稱矩陣的特征值和特征向量

2.2 主要過程

   特征臉EigenFace的思想是把人臉從像素空間變換到另一個空間，在另一個空間中做相似性計算，EigenFace選擇的空間變換方法是PCA，就是大名鼎鼎的主成分分析。EigenFace方法利用PCA得到人臉分布的主成分，具體實現是對訓練集中的所有人臉圖像的協方差矩陣進行求特征值，特征值對應的特征向量就是所謂的“特征臉”，每個特征向量描述人臉的一種變化或者特征，所以每個人臉都可以表示為這些特征臉的線性組合。下面結合以AT&T人臉庫（40個人每個人包含10個表情臉圖像，共400個臉部圖像，每個圖像分辨率為92x112），取其中399個人臉為樣本庫，最后1個為待識別人臉，給出基於Eigen特征臉的人臉識別實現過程：

（1）將訓練集中的每一個人臉圖像數據都拉長成一行，並將他們組合在一起形成一個大矩陣A，則A的大小為399x10304，即399行10304列。

（2）將399個人臉每個人臉對應的維度數據相加，然后求平均值，得到平均值向量Mean_1x10304，將矩陣A的每一行都減去平均值向量得到差值矩陣B。

（3）計算協方差矩陣C=B*B^T，C的維度是399x399，再對C求特征值λ_i，及特征向量e_i，0<=i<399。

（4）上一步驟中其實並不是真正的人臉取樣集協方差矩陣，因為人臉取樣的維度是10304，而協方差矩陣反應的是各個維度之前的相關性，所以人臉取樣集真正的協方差矩陣是C'=C^T=B^T*B，如果v_i是C'的第i個特征向量，可以證明λ_i同樣是C'的特征值，且v_i=B^T*e_i（v_i是10304行列向量），證明如下：

C*e_i=λ_i*e_i  =>  B*B^T*e_i=λ_i*e_i  =>  B^T*B*B^T*e_i=λ_i*B^T*e_i  =>  C'*v_i=λ_i*v_i 

特征向量v_i即為“特征臉”，所有特征向量組成特征向量矩陣V_10304*399，則對於任意人臉向量α，將它與特征向量矩陣V相乘，將得到向量α在各個特征向量的投影，即α*V所得向量的每一個元素為α在對應“特征臉”的投影，在進行識別時，先求得待識別人臉向量在“特征臉”的投影向量，之后和每個樣本臉的投影向量進行相似度比較，相似度最低者為最佳匹配。

2.3 核心源碼

代碼取自Opencv2.4.9

void Eigenfaces::train(InputArrayOfArrays _src, InputArray _local_labels) {
    if(_src.total() == 0) {
        string error_message = format("Empty training data was given. You'll need more than one sample to learn a model.");
        CV_Error(CV_StsBadArg, error_message);
    } else if(_local_labels.getMat().type() != CV_32SC1) {
        string error_message = format("Labels must be given as integer (CV_32SC1). Expected %d, but was %d.", CV_32SC1, _local_labels.type());
        CV_Error(CV_StsBadArg, error_message);
    }
    // make sure data has correct size
    if(_src.total() > 1) {
        for(int i = 1; i < static_cast<int>(_src.total()); i++) {
            if(_src.getMat(i-1).total() != _src.getMat(i).total()) {
                string error_message = format("In the Eigenfaces method all input samples (training images) must be of equal size! Expected %d pixels, but was %d pixels.", _src.getMat(i-1).total(), _src.getMat(i).total());
                CV_Error(CV_StsUnsupportedFormat, error_message);
            }
        }
    }
    // get labels
    Mat labels = _local_labels.getMat();
    // observations in row
    Mat data = asRowMatrix(_src, CV_64FC1);

    // number of samples
   int n = data.rows;
    // assert there are as much samples as labels
    if(static_cast<int>(labels.total()) != n) {
        string error_message = format("The number of samples (src) must equal the number of labels (labels)! len(src)=%d, len(labels)=%d.", n, labels.total());
        CV_Error(CV_StsBadArg, error_message);
    }
    // clear existing model data
    _labels.release();
    _projections.clear();
    // clip number of components to be valid
    if((_num_components <= 0) || (_num_components > n))
        _num_components = n;

    // perform the PCA
    PCA pca(data, Mat(), CV_PCA_DATA_AS_ROW, _num_components);
    // copy the PCA results
    _mean = pca.mean.reshape(1,1); // store the mean vector
    _eigenvalues = pca.eigenvalues.clone(); // eigenvalues by row
    transpose(pca.eigenvectors, _eigenvectors); // eigenvectors by column
    // store labels for prediction
    _labels = labels.clone();
    // save projections
    for(int sampleIdx = 0; sampleIdx < data.rows; sampleIdx++) {
        Mat p = subspaceProject(_eigenvectors, _mean, data.row(sampleIdx));
        _projections.push_back(p);
    }
}

 38行的PCA類中實現了求樣本矩陣的協方差矩陣、求協方差矩陣特征向量等核心功能，47行_mean為人臉平均值向量，該行其實是求每一個人臉向量減去平均值向量在“特征臉”集上的投影向量。

PCA& PCA::operator()(InputArray _data, InputArray __mean, int flags, int maxComponents)
{
    Mat data = _data.getMat(), _mean = __mean.getMat();
    int covar_flags = CV_COVAR_SCALE;
    int i, len, in_count;
    Size mean_sz;

    CV_Assert( data.channels() == 1 );
    if( flags & CV_PCA_DATA_AS_COL )
    {
        len = data.rows;
        in_count = data.cols;
        covar_flags |= CV_COVAR_COLS;
        mean_sz = Size(1, len);
    }
    else
    {
        len = data.cols;
        in_count = data.rows;
        covar_flags |= CV_COVAR_ROWS;
        mean_sz = Size(len, 1);
    }

    int count = std::min(len, in_count), out_count = count;
    if( maxComponents > 0 )
        out_count = std::min(count, maxComponents);

    // "scrambled" way to compute PCA (when cols(A)>rows(A)):
    // B = A'A; B*x=b*x; C = AA'; C*y=c*y -> AA'*y=c*y -> A'A*(A'*y)=c*(A'*y) -> c = b, x=A'*y
    if( len <= in_count )
        covar_flags |= CV_COVAR_NORMAL;

    int ctype = std::max(CV_32F, data.depth());
    mean.create( mean_sz, ctype );

    Mat covar( count, count, ctype );

    if( _mean.data )
    {
        CV_Assert( _mean.size() == mean_sz );
        _mean.convertTo(mean, ctype);
        covar_flags |= CV_COVAR_USE_AVG;
    }

    calcCovarMatrix( data, covar, mean, covar_flags, ctype );
    eigen( covar, eigenvalues, eigenvectors );

    if( !(covar_flags & CV_COVAR_NORMAL) )
    {
        // CV_PCA_DATA_AS_ROW: cols(A)>rows(A). x=A'*y -> x'=y'*A
        // CV_PCA_DATA_AS_COL: rows(A)>cols(A). x=A''*y -> x'=y'*A'
        Mat tmp_data, tmp_mean = repeat(mean, data.rows/mean.rows, data.cols/mean.cols);
        if( data.type() != ctype || tmp_mean.data == mean.data )
        {
            data.convertTo( tmp_data, ctype );
            subtract( tmp_data, tmp_mean, tmp_data );
        }
        else
        {
            subtract( data, tmp_mean, tmp_mean );
            tmp_data = tmp_mean;
        }

        Mat evects1(count, len, ctype);
        gemm( eigenvectors, tmp_data, 1, Mat(), 0, evects1,
            (flags & CV_PCA_DATA_AS_COL) ? CV_GEMM_B_T : 0);
        eigenvectors = evects1;

        // normalize eigenvectors
        for( i = 0; i < out_count; i++ )
        {
            Mat vec = eigenvectors.row(i);
            normalize(vec, vec);
        }
    }

    if( count > out_count )
    {
        // use clone() to physically copy the data and thus deallocate the original matrices
        eigenvalues = eigenvalues.rowRange(0,out_count).clone();
        eigenvectors = eigenvectors.rowRange(0,out_count).clone();
    }
    return *this;
}

45行求樣本矩陣的協方差矩陣，46行求協方差矩陣的特征值及特征向量。

void Eigenfaces::predict(InputArray _src, int &minClass, double &minDist) const {
    // get data
    Mat src = _src.getMat();
    // make sure the user is passing correct data
    if(_projections.empty()) {
        // throw error if no data (or simply return -1?)
        string error_message = "This Eigenfaces model is not computed yet. Did you call Eigenfaces::train?";
        CV_Error(CV_StsError, error_message);
    } else if(_eigenvectors.rows != static_cast<int>(src.total())) {
        // check data alignment just for clearer exception messages
        string error_message = format("Wrong input image size. Reason: Training and Test images must be of equal size! Expected an image with %d elements, but got %d.", _eigenvectors.rows, src.total());
        CV_Error(CV_StsBadArg, error_message);
    }
    // project into PCA subspace
    Mat q = subspaceProject(_eigenvectors, _mean, src.reshape(1,1));
    minDist = DBL_MAX;
    minClass = -1;
    for(size_t sampleIdx = 0; sampleIdx < _projections.size(); sampleIdx++) {
        double dist = norm(_projections[sampleIdx], q, NORM_L2);
        if((dist < minDist) && (dist < _threshold)) {
            minDist = dist;
            minClass = _labels.at<int>((int)sampleIdx);
        }
    }
}

 15行求待識別人臉向量減去人臉平均值向量在“特征臉”集上的投影向量X，19行求X與人臉樣本投影向量的歐幾里得距離（把此距離作為人臉相似度），20~23行取最小距離為識別結果。

3 示例代碼

最后給出Eigen人臉識別的示例代碼，代碼中仍使用AT&T人臉庫，其下載地址見上一篇隨筆。

#include "opencv2/core/core.hpp"
#include "opencv2/highgui/highgui.hpp"
#include "opencv2/contrib/contrib.hpp"

#define CV_VERSION_ID       CVAUX_STR(CV_MAJOR_VERSION) CVAUX_STR(CV_MINOR_VERSION) CVAUX_STR(CV_SUBMINOR_VERSION)

#ifdef _DEBUG
#define cvLIB(name) "opencv_" name CV_VERSION_ID "d"
#else
#define cvLIB(name) "opencv_" name CV_VERSION_ID
#endif

#pragma comment( lib, cvLIB("core") )
#pragma comment( lib, cvLIB("imgproc") )
#pragma comment( lib, cvLIB("highgui") )
#pragma comment( lib, cvLIB("flann") )
#pragma comment( lib, cvLIB("features2d") )
#pragma comment( lib, cvLIB("calib3d") )
#pragma comment( lib, cvLIB("gpu") )
#pragma comment( lib, cvLIB("legacy") )
#pragma comment( lib, cvLIB("ml") )
#pragma comment( lib, cvLIB("objdetect") )
#pragma comment( lib, cvLIB("ts") )
#pragma comment( lib, cvLIB("video") )
#pragma comment( lib, cvLIB("contrib") )
#pragma comment( lib, cvLIB("nonfree") )

#include <iostream>
#include <fstream>
#include <sstream>

using namespace cv;
using namespace std;

static Mat toGrayscale(InputArray _src) {
    Mat src = _src.getMat();
    // only allow one channel
    if(src.channels() != 1) {
        CV_Error(CV_StsBadArg, "Only Matrices with one channel are supported");
    }
    // create and return normalized image
    Mat dst;
    cv::normalize(_src, dst, 0, 255, NORM_MINMAX, CV_8UC1);
    return dst;
}

static void read_csv(const string& filename, vector<Mat>& images, vector<int>& labels, char separator = ';') {
    std::ifstream file(filename.c_str(), ifstream::in);
    if (!file) {
        string error_message = "No valid input file was given, please check the given filename.";
        CV_Error(CV_StsBadArg, error_message);
    }
    string line, path, classlabel;
    while (getline(file, line)) {
        stringstream liness(line);
        getline(liness, path, separator);
        getline(liness, classlabel);
        if(!path.empty() && !classlabel.empty()) {
            images.push_back(imread(path, 0));
            labels.push_back(atoi(classlabel.c_str()));
        }
    }
}

int main(int argc, const char *argv[]) {
    // Check for valid command line arguments, print usage
    // if no arguments were given.
    if (argc != 2) {
        cout << "usage: " << argv[0] << " <csv.ext>" << endl;
        exit(1);
    }

    // Get the path to your CSV.
    string fn_csv = string(argv[1]);
    // These vectors hold the images and corresponding labels.
    vector<Mat> images;
    vector<int> labels;
    // Read in the data. This can fail if no valid
    // input filename is given.
    try {
        read_csv(fn_csv, images, labels);
    } catch (cv::Exception& e) {
        cerr << "Error opening file \"" << fn_csv << "\". Reason: " << e.msg << endl;
        // nothing more we can do
        exit(1);
    }
    // Quit if there are not enough images for this demo.
    if(images.size() <= 1) {
        string error_message = "This demo needs at least 2 images to work. Please add more images to your data set!";
        CV_Error(CV_StsError, error_message);
    }
    // Get the height from the first image. We'll need this
    // later in code to reshape the images to their original
    // size:
    int height = images[0].rows;
    // The following lines simply get the last images from
    // your dataset and remove it from the vector. This is
    // done, so that the training data (which we learn the
    // cv::FaceRecognizer on) and the test data we test
    // the model with, do not overlap.
    Mat testSample = images[images.size() - 1];
    int testLabel = labels[labels.size() - 1];
    images.pop_back();
    labels.pop_back();
    // The following lines create an Eigenfaces model for
    // face recognition and train it with the images and
    // labels read from the given CSV file.
    // This here is a full PCA, if you just want to keep
    // 10 principal components (read Eigenfaces), then call
    // the factory method like this:
    //
    //      cv::createEigenFaceRecognizer(10);
    //
    // If you want to create a FaceRecognizer with a
    // confidennce threshold, call it with:
    //
    //      cv::createEigenFaceRecognizer(10, 123.0);
    //
    Ptr<FaceRecognizer> model = createEigenFaceRecognizer();
    model->train(images, labels);
    // The following line predicts the label of a given
    // test image:
    int predictedLabel = model->predict(testSample);
    //
    // To get the confidence of a prediction call the model with:
    //
    //      int predictedLabel = -1;
    //      double confidence = 0.0;
    //      model->predict(testSample, predictedLabel, confidence);
    //
    string result_message = format("Predicted class = %d / Actual class = %d.", predictedLabel, testLabel);
    cout << result_message << endl;
    // Sometimes you'll need to get/set internal model data,
    // which isn't exposed by the public cv::FaceRecognizer.
    // Since each cv::FaceRecognizer is derived from a
    // cv::Algorithm, you can query the data.
    //
    // First we'll use it to set the threshold of the FaceRecognizer
    // to 0.0 without retraining the model. This can be useful if
    // you are evaluating the model:
    //
    model->set("threshold", 0.0);
    // Now the threshold of this model is set to 0.0. A prediction
    // now returns -1, as it's impossible to have a distance below
    // it
    predictedLabel = model->predict(testSample);
    cout << "Predicted class = " << predictedLabel << endl;
    // Here is how to get the eigenvalues of this Eigenfaces model:
    Mat eigenvalues = model->getMat("eigenvalues");
    // And we can do the same to display the Eigenvectors (read Eigenfaces):
    Mat W = model->getMat("eigenvectors");
    // From this we will display the (at most) first 10 Eigenfaces:
    for (int i = 0; i < min(10, W.cols); i++) {
        string msg = format("Eigenvalue #%d = %.5f", i, eigenvalues.at<double>(i));
        cout << msg << endl;
        // get eigenvector #i
        Mat ev = W.col(i).clone();
        // Reshape to original size & normalize to [0...255] for imshow.
        Mat grayscale = toGrayscale(ev.reshape(1, height));
        // Show the image & apply a Jet colormap for better sensing.
        Mat cgrayscale;
        applyColorMap(grayscale, cgrayscale, COLORMAP_JET);
        imshow(format("%d", i), cgrayscale);
    }
    waitKey(0);

    return 0;
}

程序運行結果及用偽彩色圖像顯示的前10個特征臉，如圖所示：

本博客參考了以下資料，一並致謝！

http://www.cnblogs.com/guoming0000/archive/2012/09/27/2706019.html

http://blog.csdn.net/zouxy09/article/details/45276053

http://blog.csdn.net/zhouxuguang236/article/details/40212143

http://wenku.baidu.com/view/6023207e168884868762d644.html

《數值分析簡明教程》 王兵團 張作泉 趙平福 編著