基于PCA实现人脸识别

---

参考网址：

基于PCA的人脸识别步骤 http://blog.csdn.net/yutianzuijin/article/details/10823985

PCA人脸识别学习及C语言实现 http://blog.csdn.net/jinshengtao/article/details/18599165

PCA(主分量分析)在人脸识别中的应用

---

主题摘要

人脸识别主要方法：

  .Eigenfaces，PCA（Principal Component Analysis），Turk and Pentland，1991  .Fisherfaces，LDA（Linear Discriminant Analysis），Belhumeur, Hespanha and Kriegman，1997  .LBPH，Local Binary Pattern Histograms，Ahonen, Hadid and Pietikäinen，2004

在这篇博客中，我们主要是讲PCA的人脸识别如何实现，关于其他方法我们之后再讲。

人脸识别是一个有监督学习过程，首先利用训练集构造一个人脸模型，然后将测试集与训练集进行匹配，找到与之对应的训练集头像。最容易的方式是直接利用欧式距离计算测试集的每一幅图像与训练集的每一幅图像的距离，然后选择距离最近的图像作为识别的结果。这种直接计算距离的方式直观，但是有一个非常大的缺陷—计算量太大。如果每幅图像大小为100*100，训练集大小1000，则识别测试集中的一幅图像就需要1000*100*100的计算量，当测试集很大时，识别速度非常缓慢。

解决上述问题的一个途径是对图像进行降维，通过只保留某些关键像素可以使识别速度大大提升。降维的一个方法即是PCA（主成分分析），在此我们介绍通过PCA进行人脸识别的步骤。

一、读取训练集图像数据

读取测试集目录下指定个数的图像，然后将其保存在一个二维数组中。如果图像个数为m，图像长宽为i、j，则我们创建一个二维数组A[m][i*j=n]用来保存图像数据。数组的每一行表示一个图像的所有像素信息，每一列表示一个随机变量，也即不同图像同一位置的像素信息，降维也即用更少的列来代表图像。

二、每列减去均值

将步骤一的每列减去该列的均值，这样每列的数据均值为0。在利用matlab的函数princomp执行PCA的过程中，princomp会首先将每一列减去该列均值，不用我们自己执行。

三、计算协方差矩阵

协方差矩阵表示不同随机变量之间的相互关系，图像中也即求任意两个像素之间的关系。如果两个随机变量的协方差为正或为负，表明两个变量之间具有相关性，如果为零表示两个变量不相关。通过计算协方差矩阵，我们就可以获得不同像素之间的关系。针对人脸识别，计算的协方差矩阵大小为n*n，其中n表示图像的像素点个数。

四、计算协方差矩阵的特征值和特征向量

由于协方差矩阵是实对称阵，所以可以求得其所有的特征值和特征向量，其共有n个特征值和特征向量。

五、选择主成份

所谓主成分即是具有最大特征值的特征向量，所以我们需要将特征向量按照特征值由大到小排序，然后根据精度要求选择不同数量的特征向量，例如我们选择了前p个特征向量，通常p远小于n（在我们的人脸识别实验中，为了达到95%的精度，p只有72，而n为120*140=17040）。

六、将训练集进行降维

此步骤将原始的训练集进行降维变换，原始的图像数据是m*n的矩阵，只包含主成分的特征向量构成一个n*p的矩阵（每一列都是一个特征向量）。将两个矩阵相乘，我们即可获得降维之后的图像矩阵m*p，这个矩阵远小于原始的图像数据。

七、将测试集进行降维

同步骤6相似，读取所有的测试集图像，然后对其也进行降维操作。如果测试集有M幅图像，则降维后的矩阵为M*p。

八、人脸识别

该步骤为人脸识别的最后一步，用来对测试集进行识别，并计算识别准确率。该步骤有一个限制，测试集中的头像必须包含在训练集中，否则得出的结果将没有意义（这也就是代码一开始要求训练集大于测试集的目的）。识别的方法和最初的图像匹配方法类似：将测试集中的每一幅降维图像与降维的训练集进行匹配，然后将其分类到距离最小的训练集头像中，如果两个头像表示一个人，表示识别成功，否则表示识别失败。与原始的匹配相比，由于对图像进行了降维，所以匹配速度大大提升，在我们的实验中速度提升了200以上（120*142/72）。
如果想了解PCA的原理，请参考：PCA的数学原理

人脸识别匹配结果如图所示：(注意需要对样本进行归一化处理)

总结

该算法有以下优缺点
1）对光线的影响不是很大。与实验结果吻合，在实验时，曾经在关闭所有灯的情况下测试也可以识别出来。
2）受姿势的影响较大，只能达到50%（该结果是根据一个论文上得来的），与我们的实验结果相近。
3）论文上说 受离镜头的远近影响更大，只能识别1/3，但是我们通过一个简单的预处理步骤，把所有的人脸放大到相同的大小，这样的识别率便会比较高。
缺点么，好象识别率需要进一步验证，用于实际中不是很稳定，（当然，使用OrL库来识别时，识别率可以达到95%左右）。
现在进一步需要完善的是，
1）针对于目前应用于实际中的识别率较低，原因可能是出于检测时不太精确，虽然也采取了一种很粗糙的方法（把框住的人脸矩形往里面缩小，使得到的人脸图片近似于ORL库中的图片），相信该方法可以在一定程度上提高识别率。
2）目前所用的特征系数，是把每个人的每张图片的系数分解出来以后，在作一个平均，这种方法虽好，但是不能表示各个不同表情，各个姿势的人，下一步的工作是每个人的每幅图片都保存一个分解系数，识别时依次检索。这样待识别的人脸有可能不近似于这个人的平均脸，但极相似于该人的某一张图片，这样检索效果可能有部分提高），另外此种方法也有弊，待识别的人脸可能近似于另外一个人的某张图片。

下面要写一个关于这方面的论文：
提纲大概如下：
1）PCA算法把高维的信息转化为低维。（对于多个人，维数是否真的会降低？例如48*48的是48^2*48^2维，而对于50个人则是50*50维，维数确实有很大的降低，但是对于2304个人呢？则维数马上又上去了！！！）
2）类内和类间 的差别。在训练集中每个人可能有多个照片，这时如何处理便是一个很大的问题，一种方法就是把求出来的分解系数进行平均，另外一种方法就像前面所说的需要完善的第二个方面。
3)人脸训练出来的数据要进行归一化，在识别之前也要进行归一化。
4）人脸识别的距离判断是欧氏距离，需要进一步的完善。如果有可能的话可以试一试其他方法，诸如SVM(支持向量机）贝叶斯分类器等等。
5）人脸识别的来源是人脸检测的延续，所以需要进一步研究人脸检测的东西，使框下来的人脸大小适中，并且包含主要信息。

PCA代码：

#include <iostream>#include "opencv2/opencv.hpp"#include "Eigen/Dense"#include <fstream>#include "dirent.h"#include "time.h"using namespace cv;using namespace std;using namespace Eigen;#define INF 9999999MatrixXd getFileData2Matrix(string fileName, int featureNum, int SampleNum);MatrixXd getImageData2Matrix(string fileName, int featureNum, int SampleNum, vector<string> &ImageName);MatrixXd PCA_Dimension(MatrixXd InputTrain,int featureNUm, int SampleNum);MatrixXd normlizationMAX_MIN(MatrixXd input, int featureNum, int SampleNum);int rows, cols;vector<string> TrainImageName;vector<string> TestImageName;int main(){    // 输入 ./dataset/dataset1/Train2 ./dataset/dataset1/Test2 360 12 4    string TrainDirName, TestDirName;    int featureNum, TrainNum, TestNum;    cout<<"请输入训练样本目录名称， 测试样本目录名称，样本图片像素 ,训练样本个数和测试样本个数"<<endl;    cin>>TrainDirName>>TestDirName>>featureNum>>TrainNum>>TestNum;    //获取Train和Test数据集    MatrixXd InputTrain = getImageData2Matrix(TrainDirName,featureNum, TrainNum, TrainImageName);    MatrixXd InputTest = getImageData2Matrix(TestDirName, featureNum, TestNum, TestImageName);    for(int i=0;i<(int)TrainImageName.size();i++) cout<<"训练集图片："<<TrainImageName[i]<<endl;    for(int i=0;i<(int)TestImageName.size();i++) cout<<"测试集图片："<<TestImageName[i]<<endl;    /*    cout<<"为新的文件数据集输入featureNum 和 TrainNum"<<endl;    cin>>featureNum>>TrainNum;    MatrixXd FileMatrix = getFileData2Matrix("input1.txt", featureNum, TrainNum);    cout<<"原始样本矩阵为："<<endl<<FileMatrix<<endl;    FileMatrix = normlizationMAX_MIN(FileMatrix, featureNum, TrainNum);    cout<<"最大最小归一化矩阵为："<<endl<<FileMatrix<<endl;    MatrixXd PCAMatrix = PCA_Dimension(FileMatrix,featureNum, TrainNum);    cout<<"投影矩阵为："<<endl<<PCAMatrix<<endl;    */    clock_t start = clock();    //最大最小值归一化    //InputTrain = normlizationMAX_MIN(InputTrain, featureNum, TrainNum);    //InputTest = normlizationMAX_MIN(InputTest, featureNum, TestNum);    //使用PCA获取投影矩阵    MatrixXd ProjectMatrix = PCA_Dimension(InputTrain,featureNum, TrainNum);    //计算降维的训练样本矩阵和测试样本矩阵    MatrixXd LowDimensionTrain = InputTrain * ProjectMatrix;    MatrixXd LowDimensionTest = InputTest * ProjectMatrix;    cout<<"降维之后的训练样本矩阵大小为："<<LowDimensionTrain.rows()<<"x"<<LowDimensionTrain.cols()<<endl;    //比较训练集数据与测试集数据相似度    for(int i=0;i<TestNum;i++)    {        double dist = INF;        int index = 0;        for(int j=0;j<TrainNum;j++)        {            double dis = 0;            for(int t=0;t<LowDimensionTrain.cols();t++)            {                dis += (  LowDimensionTest(i,t) - LowDimensionTrain(j,t)  )  * (  LowDimensionTest(i,t) - LowDimensionTrain(j,t)  );            }            if(dis < dist)            {                dist = dis;                index = j;            }        }        cout<<"测试集的"<<TestImageName[i]<<"图片与训练集的"<<TrainImageName[index]<<"图片匹配度最高!!!!"<<endl;    }    clock_t end = clock();    printf("程序运行时间为:%fs\n", (float)(end - start) / CLOCKS_PER_SEC);    return 0;}MatrixXd normlizationMAX_MIN(MatrixXd input, int featureNum, int SampleNum){    MatrixXd inputMax(1,featureNum);    MatrixXd inputMin(1,featureNum);    //初始化    for(int i=0;i<featureNum;i++)    {        inputMax(0,i) = 0;        inputMin(0,i) = INF;    }    //寻找最大最小值    for(int i=0;i<featureNum;i++)    {        for(int j=0;j<SampleNum;j++)        {            if( inputMax(0,i) < input(j,i))            inputMax(0,i) = input(j,i);            if( inputMin(0,i) > input(j,i))            inputMin(0,i) = input(j,i);        }    }    //归一化    for(int i=0;i<SampleNum;i++)    {        for(int j=0;j<featureNum;j++)        {            input(i,j) = ( input(i,j) - inputMin(0,j) ) / ( inputMax(0,j) - inputMin(0,j) );        }    }    return input;}MatrixXd PCA_Dimension(MatrixXd InputTrain,int featureNum, int SampleNum){    MatrixXd CovMatrix(featureNum, featureNum);      MatrixXd MeanMatrix(1,featureNum);    double maxRate;    cout<<"请输入理想的贡献率"<<endl;    cin>>maxRate;    for(int j=0;j<featureNum;j++)    {        MeanMatrix(0,j) = 0;        for(int i=0;i<SampleNum;i++) MeanMatrix(0,j) += InputTrain(i,j);        MeanMatrix(0,j) /= SampleNum;    }    //cout<<"均值矩阵为："<<endl<<MeanMatrix<<endl;    Mat MeanFace = Mat::zeros(rows, cols, CV_8UC1);    for(int i=0;i<rows;i++)    {        for(int j=0;j<cols;j++)        MeanFace.at<uchar>(i,j) = MeanMatrix(0, i * cols + j);    }    imwrite("MeanFace.jpg",MeanFace);    double num =0;    for(int i=0;i<featureNum;i++)    {        for(int j=0;j<featureNum;j++)        {            num = 0;            for(int t=0;t<SampleNum;t++)            {                num += (InputTrain(t,i) - MeanMatrix(i)) * (InputTrain(t,j) - MeanMatrix(j));            }            num = num / (SampleNum - 1);            CovMatrix(i,j) = num;        }    }    //cout<<"协方差矩阵为："<<endl<<CovMatrix<<endl;    //计算特征值 和 特征向量    EigenSolver<MatrixXd> es(CovMatrix);    MatrixXd V = es.pseudoEigenvectors();    MatrixXd D = es.pseudoEigenvalueMatrix();    //cout<<"特征向量V:"<<endl<<V<<endl;    //cout<<"特征值矩阵D:"<<endl<<D<<endl;    // 将特征值从小到大排列    multimap<double,int> mapEigen;    for(int i=0;i<featureNum;i++)    {        pair<double, int>p(D(i,i), i);        mapEigen.insert(p);    }    //根据特征值大小 排序 特征矩阵    MatrixXd TempMatrix(featureNum, featureNum);    multimap<double, int>::iterator it = mapEigen.end();    MatrixXd EigenValue(1, featureNum);    int col = 0;    while(it != mapEigen.begin())    {        it--;        for(int i=0;i<featureNum;i++)        {            TempMatrix(i, col) = V(i, (*it).second);        }        EigenValue(0,col) = (*it).first;        col++;    }     //cout<<"按从大到小排序的特征值："<<endl<<EigenValue<<endl;    //cout<<"根据特征值排好序的特征矩阵："<<endl<<TempMatrix<<endl;    //利用特征值计算贡献率    double rateSum1 = 0, rateSum2 = 0, rate = 0;    int k;    for(int i=0;i<featureNum;i++) rateSum1 += EigenValue(0,i);    for(int i=0;i<featureNum;i++)    {        rateSum2 += EigenValue(0,i);        rate = rateSum2 / rateSum1;        cout<<"rate = "<<rate<<endl;        if(rate >= maxRate)        {            k = i + 1;            break;        }    }    MatrixXd ProjectMatrix(featureNum, k);    for(int i=0;i<featureNum;i++)    for(int j=0;j<k;j++)        ProjectMatrix(i,j) = TempMatrix(i,j);    cout<<"前"<<k<<"个"<<"贡献率 >= "<<maxRate<<"的投影矩阵"<<endl;    return ProjectMatrix;}//读取样本数据集MatrixXd getFileData2Matrix(string fileName, int featureNum, int SampleNum){    MatrixXd InputTrain(SampleNum,featureNum);    cout<<"矩阵大小为："<<InputTrain.rows()<<"x"<<InputTrain.cols()<<endl;    ifstream in(fileName);    if( in.is_open())    {        for(int i=0;i<SampleNum;i++)        {            for(int j=0;j<featureNum;j++)            {                in>>InputTrain(i,j);            }        }        return InputTrain;    }    else     {         cout<<"打开文件失败，文件名称输入不正确或文件不存在"<<endl;         return InputTrain;    }}MatrixXd getImageData2Matrix(string fileName, int featureNum, int SampleNum, vector<string> &ImageName){    Mat pic;    MatrixXd InputTrain(SampleNum, featureNum);    cout<<"矩阵大小为"<<InputTrain.rows()<<"x"<<InputTrain.cols()<<endl;    DIR* dir = opendir((char*)fileName.c_str());//打开指定目录      dirent* p = NULL;//定义遍历指针      int index = 0;    while((p = readdir(dir)) != NULL)//开始逐个遍历      {          //这里需要注意，linux平台下一个目录中有"."和".."隐藏文件，需要过滤掉          if(p->d_name[0] != '.')//d_name是一个char数组，存放当前遍历到的文件名          {              string name = fileName +"/"+ string(p->d_name);              ImageName.push_back( string(p->d_name) );            //cout<<name<<endl;            pic = imread(name);            cvtColor(pic,pic, CV_RGB2GRAY);            if(!pic.isContinuous())            {                cout<<"样本图片数据不连续"<<endl;                return InputTrain;            }            rows = pic.rows, cols = pic.cols;            cout<<"图片大小为"<<rows<<"x"<<cols<<endl;            for(int i = 0;i<rows;i++)            {                for(int j = 0;j<cols;j++)                {                    InputTrain(index, i*cols+j) = (double)pic.at<uchar>(i,j);                    //cout<<i<<" "<<j<<endl;                }            }            index++;            //cout<<"index = "<<index<<endl;        }      }      closedir(dir);//关闭指定目录      return InputTrain;}

如果输入的图片的尺寸太大，那么应在保证尽量不改变图片信息的前提下，适当地缩小图片尺寸。
ChangeSize:

#include <iostream>#include "opencv2/opencv.hpp"#include "dirent.h"using namespace std;using namespace cv;void ChangeSize(string sourceDir, string saveDir){    DIR* dir = opendir((char*)sourceDir.c_str());//打开指定目录      dirent* p = NULL;//定义遍历指针      Mat pic;    while((p = readdir(dir)) != NULL)//开始逐个遍历      {          //这里需要注意，linux平台下一个目录中有"."和".."隐藏文件，需要过滤掉          if(p->d_name[0] != '.')//d_name是一个char数组，存放当前遍历到的文件名          {              string name = sourceDir +"/"+ string(p->d_name);              //cout<<name<<endl;            pic = imread(name);            cvtColor(pic,pic, CV_RGB2GRAY);            resize(pic, pic, Size(28, 23));            imshow("Pic",pic);            string str = saveDir + "/" +string(p->d_name);            imwrite(str, pic);            waitKey(50);        }      }      closedir(dir);//关闭指定目录      return;}int main(){    string sourceDir, saveDir;    cout<<"请输入图片来源目录名称 和 保存目录名称 "<<endl;    cin>>sourceDir>>saveDir;    ChangeSize(sourceDir, saveDir);}