（学习笔记三）——人脸识别

原文地址：http://blog.csdn.net/longxiaoshi/article/details/7411333

人脸识别要牵涉到一些数学计算和一些算法的理解，虽然这些算法和计算opencv已经帮我们完成，但我们还是要对其有一定的了解，才能进行人脸识别的实践，毕竟基础不牢，上层建筑也不稳。

要理解如何进行人脸识别，首先一定要理解主成分分析算法，即PCA，使用这种算法的原因是因为，一般图像数据量太大，而且其中的大部分数据点对我们进行人脸识别没有太大的帮助，因此为了减少数据量，采用了主成分分析法将数据进行压缩（算法中称为投影），然后用压缩后的图像进行图像识别的进一步应用。PCA的原理在前面的一篇文章中已经讲到过，大家可以参考： PCA原理。肯能看完这个文章，大家开始对PCA有一定的理解，但总还是感觉有点模糊，因为那篇文章我只是转载，并没有进行太多自己的想法的标注，因此有些混乱，那么大家可以再继续读下面这篇文章：
PCA原理详解。
在前面的PCA原理介绍中曾提到了两个函数，cvCalcPCA，cvProjectPCA。我在实际使用中参考了另外一篇文章，用的是另外两个函数，这个后面会讲到。这个不用太纠结，前面的文章可以只当学习原理的参考，等原理懂了之后看后面将要介绍的两个函数也会很简单。其中，有一点要说明，cvCalcPCA是新的PCA处理函数，而后面将要讲到的cvCalcEigenObjects是老的PCA处理函数。
经过上面的的两片文章相信大家应该对PCA的原理有所了解，下面开始讲怎样实现用PCA算法进行人脸识别。这个可以参考这篇文章：
http://www.cognotics.com/opencv/servo_2007_series/part_5/index.html
这篇是英文文章，希望大家多读英文文章，因为我们在学习过程中如果遇到国内教材较少的情况，一般都是要参考国外的文章的，并且擅长阅读英文文章也为我们打开了通往另外一个更为博大的知识库的大门。读完这篇文章一定可以自己做出人脸识别程序的，因为这篇文章介绍的已经十分详细。
但是，在这里我还是将上面那篇文章的原理介绍一下。我下面的叙述主要是针对在一个图像库中找出我们给定的图像最接近的图像，以人脸图像为例也就是说数据库中肯定是包含我所给出的这个人的人脸，虽然数据库中的图片跟我给出的图片不是同一张图片，但必须是同一个人。
首先假如我有代号为1,2,3三个训练人的人脸图像，人脸大小调整为一致，例如文章中为92*112，在我们的分析中是以一个像素点作为一维，因此每个图像可以认为是一个1*10304的行向量，这样三个训练人的人脸图像就组成了一个3*10304的矩阵。（注意，实际中每个图像还是存在92*112的iplimage结构中，这里只是为了讲解方便而进行的假象）。根据这个矩阵，就可以计算出相应的协方差矩阵，大小为10304*10304（这个大小我并没有考究，我觉得应该是这个大小，因为可以解释通上面那篇英文文章），接着再求出这个10304*10304的矩阵的特征值和特征向量，并取前nEigens个（这个个数由我们自己确定，文章中是取nEigens为训练人个数-1）特征值（注意由于前面特征值已经是由大到小排好顺序的，因此这几个就是最大的几个）对应的特征向量组成10304*nEigens大小投影矩阵（注意这里发生了转置，有opencv函数内部完成），即所谓的特征脸（因为每一列的规模就是一张人脸图像的规模），也即子空间。然后再将每个大小为1*10304的训练人脸图像矩阵乘以这个投影矩阵，就可以得到每个图像在主成分子空间的投影，大小为1*nEigens，并用这个投影矩阵进行后续的分析。可见，经过投影后图像数据量大大减小。上面的计算投影矩阵的过程由cvCalcEigenObjects完成，而投影则由cvEigenDecomposite函数完成，具体实现见英文文章。为了方便理解这两个函数可以参见： PCA的两个主要函数，这个空间内有两篇介绍这两个函数的文章。
同样，当给出测试图像时，先将测试人脸缩放为与训练人脸相同的大小，即92*112，然后用前面的特征矩阵将其投影到子空间中，即投影为大小为1*nEigens的矩阵。接着就可以用这个子空间投影跟前面的训练图像在子空间的投影进行比较，最接近的就是目标图像。比较方法文章中也有提到，有欧式算法 Euclidean 和较新的 Mahalanobis算法。
至此，整个人脸识别原理已经讲完，虽然有点复杂，但花点时间也不难理解。

下面贴出文章中给出的程序，我进行了一些注释，使用方法程序头部也有介绍。

// eigenface.c, by Robin Hewitt, 2007
//
// Example program showing how to implement eigenface with OpenCV

// Usage:
//
// First, you need some face images. I used the ORL face database.
// You can download it for free at
//    www.cl.cam.ac.uk/research/dtg/attarchive/facedatabase.html
//
// List the training and test face images you want to use in the
// input files train.txt and test.txt. (Example input files are provided
// in the download.) To use these input files exactly as provided, unzip
// the ORL face database, and place train.txt, test.txt, and eigenface.exe
// at the root of the unzipped database.
//
// To run the learning phase of eigenface, enter
//    eigenface train
// at the command prompt. To run the recognition phase, enter
//    eigenface test

#include <stdio.h>
#include <string.h>
#include "cv.h"
#include "cvaux.h"
#include "highgui.h"


定义几个重要的全局变量
IplImage ** faceImgArr        = 0; // 指向训练人脸和测试人脸的指针（在学习和识别阶段指向不同）
CvMat    *  personNumTruthMat = 0; // 人脸图像的ID号
int nTrainFaces               = 0; // 训练图像的数目
int nEigens                   = 0; // 自己取的主要特征值数目
IplImage * pAvgTrainImg       = 0; // 训练人脸数据的平均值
IplImage ** eigenVectArr      = 0; // 投影矩阵，也即主特征向量
CvMat * eigenValMat           = 0; // 特征值
CvMat * projectedTrainFaceMat = 0; // 训练图像的投影


 函数原型
void learn();
void recognize();
void doPCA();
void storeTrainingData();
int  loadTrainingData(CvMat ** pTrainPersonNumMat);
int  findNearestNeighbor(float * projectedTestFace);
int  loadFaceImgArray(char * filename);
void printUsage();



//主函数，主要包括学习和识别两个阶段，需要运行两次，通过命令行传入的参数区分
void main( int argc, char** argv )
{
	// validate that an input was specified
	if( argc != 2 )
	{
		printUsage();
		return;
	}
	//通过判断命令行参数分别执行学习和识别代码
	if( !strcmp(argv[1], "train") ) learn();
	else if( !strcmp(argv[1], "test") ) recognize();
	else
	{
		printf("Unknown command: %s\n", argv[1]);
		printUsage();
	}
}


//学习阶段代码
void learn()
{
	int i, offset;

	//加载训练图像集
	nTrainFaces = loadFaceImgArray("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 );
	offset = projectedTrainFaceMat->step / sizeof(float);
	for(i=0; i<nTrainFaces; i++)
	{
		//int offset = i * nEigens;
		cvEigenDecomposite(
			faceImgArr[i],
			nEigens,
			eigenVectArr,
			0, 0,
			pAvgTrainImg,
			//projectedTrainFaceMat->data.fl + i*nEigens);
			projectedTrainFaceMat->data.fl + i*offset);
	}

	//将训练阶段得到的特征值，投影矩阵等数据存为.xml文件，以备测试时使用
	storeTrainingData();
}


//识别阶段代码
void recognize()
{
	int i, nTestFaces  = 0;         // 测试人脸数
	CvMat * trainPersonNumMat = 0;  // 训练阶段的人脸数
	float * projectedTestFace = 0;

	// 加载测试图像，并返回测试人脸数
	nTestFaces = loadFaceImgArray("test.txt");
	printf("%d test faces loaded\n", nTestFaces);

	// 加载保存在.xml文件中的训练结果
	if( !loadTrainingData( &trainPersonNumMat ) ) return;

	// 
	projectedTestFace = (float *)cvAlloc( nEigens*sizeof(float) );
	for(i=0; i<nTestFaces; i++)
	{
		int iNearest, nearest, truth;

		//将测试图像投影到子空间中
		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;

	
	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);
	}

	
	cvReleaseFileStorage( &fileStorage );

	return 1;
}



//存储训练结果
void storeTrainingData()
{
	CvFileStorage * fileStorage;
	int i;

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

	//存储特征值，投影矩阵，平均矩阵等训练结果
	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));
	}


	cvReleaseFileStorage( &fileStorage );
}



//寻找最接近的图像
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 / eigenValMat->data.fl[i];  // Mahalanobis算法计算的距离
		//	distSq += d_i*d_i; // Euclidean算法计算的距离
		}

		if(distSq < leastDistSq)
		{
			leastDistSq = distSq;
			iNearest = iTrain;
		}
	}

	return iNearest;
}



//主成分分析
void doPCA()
{
	int i;
	CvTermCriteria calcLimit;
	CvSize faceImgSize;

	// 自己设置主特征值个数
	nEigens = nTrainFaces-1;

	//分配特征向量存储空间
	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);

	//分配主特征值存储空间
	eigenValMat = cvCreateMat( 1, nEigens, CV_32FC1 );

	// 分配平均图像存储空间
	pAvgTrainImg = cvCreateImage(faceImgSize, IPL_DEPTH_32F, 1);

	// 设定PCA分析结束条件
	calcLimit = cvTermCriteria( CV_TERMCRIT_ITER, nEigens, 1);

	// 计算平均图像，特征值，特征向量
	cvCalcEigenObjects(
		nTrainFaces,
		(void*)faceImgArr,
		(void*)eigenVectArr,
		CV_EIGOBJ_NO_CALLBACK,
		0,
		0,
		&calcLimit,
		pAvgTrainImg,
		eigenValMat->data.fl);

	cvNormalize(eigenValMat, eigenValMat, 1, 0, CV_L1, 0);
}



//加载txt文件的列举的图像
int loadFaceImgArray(char * filename)
{
	FILE * imgListFile = 0;
	char imgFilename[512];
	int iFace, nFaces=0;


	if( !(imgListFile = fopen(filename, "r")) )
	{
		fprintf(stderr, "Can\'t open file %s\n", filename);
		return 0;
	}

	// 统计人脸数
	while( fgets(imgFilename, 512, imgListFile) ) ++nFaces;
	rewind(imgListFile);

	// 分配人脸图像存储空间和人脸ID号存储空间
	faceImgArr        = (IplImage **)cvAlloc( nFaces*sizeof(IplImage *) );
	personNumTruthMat = cvCreateMat( 1, nFaces, CV_32SC1 );

	for(iFace=0; iFace<nFaces; iFace++)
	{
		// 从文件中读取序号和人脸名称
		fscanf(imgListFile,
			"%d %s", personNumTruthMat->data.i+iFace, imgFilename);

		// 加载人脸图像
		faceImgArr[iFace] = cvLoadImage(imgFilename, CV_LOAD_IMAGE_GRAYSCALE);

		if( !faceImgArr[iFace] )
		{
			fprintf(stderr, "Can\'t load image from %s\n", imgFilename);
			return 0;
		}
	}

	fclose(imgListFile);

	return nFaces;
}



//
void printUsage()
{
	printf("Usage: eigenface <command>\n",
		"  Valid commands are\n"
		"    train\n"
		"    test\n");
}