PCA原理及人脸识别

原文地址：  http://blog.sina.com.cn/s/blog_60a0e97e01017o5x.html

PCA原理

（原文出处：http://blog.csdn.net/longxiaoshi/article/details/7404251 ）

 

    PCA是主成分分析，主要用于数据降维，对于一系列sample的feature组成的多维向量，多维向量里的某些元素本身没有区分性，比如某个元素在所有的sample中都为1，或者与1差距不大，那么这个元素本身就没有区分性，用它做特征来区分，贡献会非常小。所以我们的目的是找那些变化大的元素，即方差大的那些维，而去除掉那些变化不大的维，从而使feature留下的都是“精品”，而且计算量也变小了。

    对于一个k维的feature来说，相当于它的每一维feature与其他维都是正交的（相当于在多维坐标系中，坐标轴都是垂直的），那么我们可以变化这些维的坐标系，从而使这个feature在某些维上方差大，而在某些维上方差很小。例如，一个45度倾斜的椭圆，在第一坐标系，如果按照x,y坐标来投影，这些点的x和y的属性很难用于区分他们，因为他们在x,y轴上坐标变化的方差都差不多，我们无法根据这个点的某个x属性来判断这个点是哪个，而如果将坐标轴旋转，以椭圆长轴为x轴，则椭圆在长轴上的分布比较长，方差大，而在短轴上的分布短，方差小，所以可以考虑只保留这些点的长轴属性，来区分椭圆上的点，这样，区分性比x,y轴的方法要好！

    所以我们的做法就是求得一个k维特征的投影矩阵，这个投影矩阵可以将feature从高维降到低维。投影矩阵也可以叫做变换矩阵。新的低维特征必须每个维都正交，特征向量都是正交的。通过求样本矩阵的协方差矩阵，然后求出协方差矩阵的特征向量，这些特征向量就可以构成这个投影矩阵了。特征向量的选择取决于协方差矩阵的特征值的大小。

 

举一个例子：

    对于一个训练集，100个sample，特征是10维，那么它可以建立一个100*10的矩阵，作为样本。求这个样本的协方差矩阵，得到一个10*10的协方差矩阵，然后求出这个协方差矩阵的特征值和特征向量，应该有10个特征值和特征向量，我们根据特征值的大小，取前四个特征值所对应的特征向量，构成一个10*4的矩阵，这个矩阵就是我们要求的特征矩阵，100*10的样本矩阵乘以这个10*4的特征矩阵，就得到了一个100*4的新的降维之后的样本矩阵，每个sample的维数下降了。

    当给定一个测试的特征集之后，比如1*10维的特征，乘以上面得到的10*4的特征矩阵，便可以得到一个1*4的特征，用这个特征去分类。

    所以做PCA实际上是求得这个投影矩阵，用高维的特征乘以这个投影矩阵，便可以将高维特征的维数下降到指定的维数。

 

    在opencv里面有专门的函数，可以得到这个这个投影矩阵（特征矩阵）。
void cvCalcPCA( const CvArr* data, CvArr* avg, CvArr* eigenvalues, CvArr* eigenvectors, int flags ); 

在OPENCV中使用PCA非常简单，只要几条语句就可以了。

• 1. 初始化数据
  //每一行表示一个样本
  CvMat* pData = cvCreateMat( 总的样本数, 每个样本的维数, CV_32FC1 );
  CvMat* pMean = cvCreateMat(1, 样本的维数, CV_32FC1);
  //pEigVals中的每个数表示一个特征值
  CvMat* pEigVals = cvCreateMat(1, min(总的样本数,样本的维数), CV_32FC1);
  //每一行表示一个特征向量
  CvMat* pEigVecs = cvCreateMat( min(总的样本数,样本的维数), 样本的维数, CV_32FC1);
• 2. PCA处理，计算出平均向量pMean,特征值pEigVals和特征向量pEigVecs
  cvCalcPCA( pData, pMean, pEigVals, pEigVecs, CV_PCA_DATA_AS_ROW );
• 3. 选出前P个特征向量(主成份)，然后投影，结果保存在pResult中，pResult中包含了P个系数
  CvMat* pResult = cvCreateMat( 总的样本数, PCA变换后的样本维数(即主成份的数目), CV_32FC1 );
  cvProjectPCA( pData, pMean, pEigVecs, pResult );
• 4. 重构，结果保存在pRecon中
  CvMat* pRecon = cvCreateMat(总的样本数, 每个样本的维数, CV_32FC1 );
  cvBackProjectPCA( pResult, pMean, pEigVecs, pRecon );
• 5.重构误差的计算
      计算pRecon和pData的"差"就可以了。
      使用时如果是想用PCA判断“是非”问题，则可以先用正样本计算主成分，判断时，对需要判断得数据进行投影，然后重构，计算重构出的数据与原数据的差异，如果差异在给定范围内，可以认为“是”。
      如果相用PCA进行分类，例如，对数字进行分类，则先用所有数据（0-9的所有样本）计算主成分，然后对每一类数据进行投影，计算投影的系数，可简单得求平均，即对每一类求出平均系数。分类时，将需要分类得数据进行投影，得到系数，与先前计算出得每一类得平均系数进行比较，可判为最接近得一类。当然这只是最简单得使用方法。

 

PCA人脸识别

（原文出处：http://www.cnblogs.com/linzhao/archive/2012/01/05/2312758.html ）

    PCA(Eigenface)方法是人脸识别的主流方法之一。cvEigenDecomposite()函数作用是将人脸图像通过Eigenface变换矩阵，投射到子空间中。子空间中的人脸向量，是一个1×nEigens（nEigens由自己取得）的行向量，极大地降低了数据维度，便于下一步的聚类、识别。结合具体例子，尽量清楚地跟大家讨论这个函数的用途。

 

函数说明如下：

EigenDecomposite：Calculates all decomposition coefficients for input object 计算输入目标的所有的分解系数

void cvEigenDecomposite( IplImage* obj, int nEigObjs, void* eigInput,int ioFlags, void* userData, IplImage* avg, float* coeffs );
 

obj - Input object，输入目标
nEigObjs - Number of eigen objects，自选eigen数量
eigInput - Pointer either to the array of IplImage input objects or to the read callback function according to the value of the parameter ioFlags，eigenface向量组成的变换矩阵，数据结构属于IplImage队列，或者根据ioFlag参数值变化的回调函数。
ioFlags - Input/output flags
userData - Pointer to the structure that contains all necessary data for the callback functions.
avg - Averaged object.
coeffs - Calculated coefficients; an output parameter 这是唯一一个输出，即人脸在子空间的投影（decomposition coefficients... 数学不太好- -）The function cvEigenDecomposite calculates all decomposition coefficients for the input object using the previously calculated eigen objects basis and the averaged object. Depending on ioFlags parameter it may be used either in direct access or callback mode.

 

    在程序中，cvCalcEigenObject()只需要在训练时使用，依据训练样本，计算Eigenface的变换矩阵。如上图所示。

    而cvEigenDecomposite()，既需要产生训练样本的投影，又需要产生测试脸的投影。

cvEigenDecomposite(faceImgArr[i], nEigens, eigenVectArr, 0, 0, pAvgTrainImg, projectedTestFace);

这里是对于每张测试脸faceImgArr[i]，产生投影。训练时得到训练样本投影脸的函数调用方法如下：

for(i=0; i<nTrainFaces; i++)
{
   //int offset = i * nEigens;
   cvEigenDecomposite(faceImgArr[i], nEigens, eigenVectArr, 0, 0, pAvgTrainImg, projectedTrainFaceMat->data.fl + i*offset);
}

 

PCA人脸识别代码实现

（原文出处：http://www.cognotics.com/opencv/servo_2007_series/part_5/index.html）

// EigenFaces.cpp : 定义控制台应用程序的入口点。

#include "stdafx.h"

#include "highgui.h"

#include "cv.h"

#include "cvaux.h"

#include "stdio.h"

#include "string.h"

int nTrainFaces =0;

int nEigens =0;

IplImage** faceImgArr =0;

CvMat* personNumTruthMat=0;

IplImage* pAvgTrainImg =0;

IplImage** eigenVectArr =0;

CvMat* eigenValMat =0;

CvMat* projectedTrainFaceMat=0;

void displayMatrix(const CvMat *mat);

void learn();

void recognize();

void doPCA();

void storeTrainingData();

int loadTrainingData(CvMat** pTrainPersonNumMat);

int findNearestNeighbor(float* projectedtestFace);

int loadFaceImgArr(char* filename);

void printUsage();

int _tmain(int argc, _TCHAR* argv[])

{

  learn();

  recognize();

  system("PAUSE");

  return 0;

}

void displayMatrix( const CvMat *mat ) 

{ 

  size_t i, j; 

  for ( i = 0; i != mat->rows; i++ ) 

  { 

    for ( j = 0; j != mat->cols; j++ ) 

    { 

      printf( "%.4lf ",cvmGet( mat, i, j ) ); 

    } 

    printf( "\n"); 

  } 

  printf( "\n" ); 

}

void learn()

{

  int i;

  nTrainFaces=loadFaceImgArr("train.txt");

  if(nTrainFaces<2)

  {

    fprintf(stderr,"Need 2 or more training faces\n","Input file contains only %d\n", nTrainFaces);

    return;

  }

  doPCA();

  projectedTrainFaceMat=cvCreateMat(nTrainFaces,nEigens,CV_32FC1);

  for (i=0;i<nTrainFaces;i++)

  {

    cvEigenDecomposite(faceImgArr[i],nEigens,eigenVectArr,0,0,

    pAvgTrainImg,projectedTrainFaceMat->data.fl+i*nEigens);

  }

  storeTrainingData();

}

int loadFaceImgArr(char *filename)

{

  FILE* imgListFile=0;

  char imgFilename[512];

  int iFace, nFaces=0;

  //open the input file

  imgListFile=fopen(filename,"r");

  while(fgets(imgFilename,512,imgListFile)) ++nFaces;

  //fgets: read a line of file's strings to store in imgFilename

  //function fgets stops if meets a line break or read number of parameter 2

  rewind(imgListFile);//get the pointer to the start position of a file

  faceImgArr=(IplImage**)cvAlloc(nFaces*sizeof(IplImage*));

  personNumTruthMat=cvCreateMat(1,nFaces,CV_32SC1);

  //store the face images in an array

  for (iFace=0;iFace<nFaces;iFace++)

  {

    //read person number and name of image file

    fscanf(imgListFile,"%d %s",personNumTruthMat->data.i+iFace,imgFilename);

    //fscanf: Reads data from the stream and stores them according to the parameter format 

    //into the locations pointed by the additional arguments. 

    //load the face image

    faceImgArr[iFace]=cvLoadImage(imgFilename,CV_LOAD_IMAGE_GRAYSCALE);

  }

  fclose(imgListFile);

  return nFaces;

}

void doPCA()

{

  int i;

  CvTermCriteria calcLimit;

  CvSize faceImgSize;

  //set the number of eigenvalues to use

  nEigens=nTrainFaces-1;

  //allocate the eigenvector image

  faceImgSize.width=faceImgArr[0]->width;

  faceImgSize.height=faceImgArr[0]->height;

  eigenVectArr=(IplImage**)cvAlloc(sizeof(IplImage*)*nEigens);

  for (i=0;i<nEigens;i++)

    eigenVectArr[i]=cvCreateImage(faceImgSize,IPL_DEPTH_32F,1);

  //allocate the eigenvalue array

  eigenValMat=cvCreateMat(1,nEigens,CV_32FC1);

  //allocate the averaged image

  pAvgTrainImg=cvCreateImage(faceImgSize,IPL_DEPTH_32F,1);

  //set the PCA termination criterion

  calcLimit=cvTermCriteria(CV_TERMCRIT_ITER,nEigens,1);

  cvCalcEigenObjects(nTrainFaces, (void*)faceImgArr, (void*)eigenVectArr,

CV_EIGOBJ_NO_CALLBACK, 0, 0, &calcLimit, pAvgTrainImg, eigenValMat->data.fl);

  printf("eigenValMat:\n");

  displayMatrix(eigenValMat);

}

void storeTrainingData()

{

  CvFileStorage* fileStorage;

  int i;

  //create a file-storage interface

  fileStorage=cvOpenFileStorage("facedata.xml",0,CV_STORAGE_WRITE);

  //store all the data

  cvWriteInt(fileStorage,"nEigens",nEigens);

  cvWriteInt(fileStorage,"nTrainFaces",nTrainFaces);

  cvWrite(fileStorage,"trainPersonNumMat",personNumTruthMat,cvAttrList(0,0));

  cvWrite(fileStorage,"eigenValMat",eigenValMat,cvAttrList(0,0));

  cvWrite(fileStorage,"projectedTrainFaceMat",projectedTrainFaceMat,cvAttrList(0,0));

  cvWrite(fileStorage,"avgTrainImg",pAvgTrainImg,cvAttrList(0,0));

  for (i=0;i<nEigens;i++)

  {

    char varname[200];

    sprintf(varname,"eigenVect_%d",i);

    cvWrite(fileStorage,varname,eigenVectArr[i],cvAttrList(0,0));

  }

  //release the file-storage interface

  cvReleaseFileStorage(&fileStorage);

}

void recognize()

{

  int i, nTestFaces=0;

  CvMat* trainPersonNumMat=0;

  float* projectedTestFace=0;

  //load test images and ground truth for person number

  nTestFaces=loadFaceImgArr("test.txt");

  printf("%d test faces loaded\n",nTestFaces);

  //load the saved training data

  if(!loadTrainingData(&trainPersonNumMat)) return;

  //project the test images onto the PCA subspace

  projectedTestFace=(float*)cvAlloc(nEigens*sizeof(float));

  for (i=0;i<nTestFaces;i++)

  {

    int iNearest, nearest, truth;

    ///project the test image onto the PCA subspace

    cvEigenDecomposite(faceImgArr[i],nEigens,eigenVectArr, 0, 0, pAvgTrainImg, projectedTestFace);

    iNearest=findNearestNeighbor(projectedTestFace);

    truth=personNumTruthMat->data.i[i];

    nearest=trainPersonNumMat->data.i[iNearest];

    printf("nearest=%d, truth=%d\n",nearest,truth);

  }

}

int loadTrainingData(CvMat** pTrainPersonNumMat)

{

  CvFileStorage* fileStorage;

  int i;

  //create a file-storage interface

  fileStorage=cvOpenFileStorage("facedata.xml",0,CV_STORAGE_READ);

  if(!fileStorage)

  {

    fprintf(stderr,"can't open facedata.xml\n");

    return 0;

  }

  nEigens=cvReadIntByName(fileStorage,0,"nEigens",0);

  nTrainFaces=cvReadIntByName(fileStorage,0,"nTrainFaces",0);

  *pTrainPersonNumMat=(CvMat*)cvReadByName(fileStorage,0,"trainPersonNumMat",0);

  eigenValMat=(CvMat*)cvReadByName(fileStorage,0,"eigenValMat",0);

  projectedTrainFaceMat=(CvMat*)cvReadByName(fileStorage,0,"projectedTrainFaceMat",0);

  pAvgTrainImg=(IplImage*)cvReadByName(fileStorage,0,"avgTrainImg",0);

  eigenVectArr=(IplImage**)cvAlloc(nTrainFaces*sizeof(IplImage*));

  for (i=0;i<nEigens;i++)

  {

    char varname[200];

    sprintf(varname,"eigenVect_%d",i);

    eigenVectArr[i]=(IplImage*)cvReadByName(fileStorage,0,varname,0);

  }

  //release the file-storage interface

  cvReleaseFileStorage(&fileStorage);

  return 1;

}

int findNearestNeighbor(float* projectedTestFace)

{

  double leastDistSq=DBL_MAX;

  int i, iTrain, iNearest=0;

  for (iTrain=0;iTrain<nTrainFaces;iTrain++)

  {

    double distSq=0;

    for (i=0;i<nEigens;i++)

    {

      float d_i=projectedTestFace[i]- projectedTrainFaceMat->data.fl[iTrain*nEigens+i];

      distSq+=d_i*d_i;

    }

    if (distSq<leastDistSq)

    {

      leastDistSq=distSq;

      iNearest=iTrain;

    }

  }

  return iNearest;

}

void printUsage()

{

  printf("Usage:EigenFaces <command>\n","Valid commands are\n","train\n","test\n");

}