• 为了保证你在浏览本网站时有着更好的体验,建议使用类似Chrome、Firefox之类的浏览器~~
    • 如果你喜欢本站的内容何不Ctrl+D收藏一下呢,与大家一起分享各种编程知识~
    • 本网站研究机器学习、计算机视觉、模式识别~当然不局限于此,生命在于折腾,何不年轻时多折腾一下

基于QR和PCA的人脸识别

cv admin 3年前 (2016-05-10) 1503次浏览 0个评论 扫描二维码

PCA(Principal Component Analysis,主成分分析)是一种很常用的根据变量协方差对数据进行降维、压缩的方法。它的精髓在于尽量用最少数量的维度,尽可能精确地描述数据。

PCA 对数据进行降维的过程可以用下面这个动图来解释(图片摘自 http://stats.stackexchange.com/a/140579/93946):

lNHqt
在上图中,一组位于直角坐标系的二维样本集,沿着斜线的方向有很强的相关性。所以如果我们将直角坐标系转换到斜向,也就是让横轴沿斜线方向,纵轴垂直于斜线方向。于是,在这个新的坐标系下,数据点在横轴上分布很分散,但是在纵轴方向比较集中。如果在误差允许范围内,我们完全可以将数据点在新纵轴上的坐标全部置为 0,只用新横轴上的坐标来表示点的位置。这样,就完成了对数据降维的过程(即将原始直角坐标系的 2 个维度减小到新坐标系的 1 个维度)。对更高维的情况,处理过程与之类似。


PCA 人脸识别

将 PCA 用于人脸识别的过程如下:

1.假设有 400 幅尺寸为 100*100 的图像,构成 10000*400 的矩阵  X=[x_1,\dots,x_n]

2.计算均值  \mu=\frac{1}{n}\sum_{j=1}^n x_j,令  H=\frac{1}{\sqrt{n-1}}[x_1-\mu,\dots,x_n-\mu]

3.根据定义,计算协方差矩阵 \Sigma=HH^T

4.计算\Sigma的特征值与特征向量,取前 h 个最大特征值所对应的特征向量,构成矩阵\Phi

5.矩阵 \Phi 可 对数据降维:\Phi^T X=Y ,Y 是 h 行 400 列的矩阵,也就是将数据从 10000 维降为 h 维。

这种做法一个明显的缺陷在于,\Sigma 的维度为 10000×10000,直接进行奇异值分解计算量非常大。利用 QR 分解,作间接的奇异值分解,可以减小计算量。


利用 QR 分解减小计算量

基于 QR 分解的 PCA 算法步骤如下:

1.已知\Sigma=HH^T,其中\Sigma为 d*d,H 为 d*n,d 代表原始数据的维数,n 代表样本数,d 远大于 n;

2.对 H 作 QR 分解,h=QR,其中 Q 为 d*t,R 为 t*n,1\leq t \leq n

3.则\Sigma=QRR^T Q^T,对R^T作奇异值分解R^T=UDV^T,其中 U 为 n*t,V 为 t*t,D=diag(\sigma_1,\dots,\sigma_t)

4.于是\Sigma=QVDU^T UDV^T Q^T=QVD^2 V^T Q^T=QV\Lambda V^T Q^T,其中\Lambda=D^2

5.由于(QV)^T (QV)=V^T Q^T QV=V^T V=I,所以 QV 可将\Sigma对角化,QV 为\Sigma的特征向量矩阵,\Lambda\Sigma的特征值矩阵;

6.选取 D 前 h 个最大对角元所对应于 V 中的 h 个列,构成 t*h 的矩阵V_h,则降维矩阵\Phi=QV_h

该过程涉及对一个很大的矩阵的 QR 分解,和对一个较小矩阵的奇异值分解。计算量与传统方法相比较的结果如下(图片摘自[1]):

computation_comparision-1024x426

进一步,进行人脸识别的过程如下:

1.假设有 c 个类别,每类包含 s 个样本,则 n=c∗s;

2.对 X 计算Y=\Phi^T X,将 Y(也称特征脸)按类别计算均值,得到 c 个长度为 h 的列向量v_1,\dots,v_c

3.对于未知类别的新样本 x,计算y=\Phi^T x,y 的长度为 h;

4.计算距离d(y,v_i),i=1,\dots,c,取距离最小的 i 作为 x 的类标号。


距离度量 d

1.欧式距离

Euclidean_Distance-300x58

2.曼哈顿距离

Manhattan_Distance-300x67

3.马氏距离

Mahalanobis_Distance

在马氏距离中,x 与 y 分布相同,且协方差矩阵为 S。加入协方差矩阵的逆矩阵的作用是,将如下图(图片部分取自 http://stats.stackexchange.com/a/62147/93946)中呈椭圆分布的数据归一化到圆形分布中,再来比较距离,可以抵消不同样本集在特征空间中的分布差异。

Mahalanobis_Distance_example


C++实现

环境要求:OpenCV(样本图像的读取),Armadillo(高性能线性代数 C++函数库),Intel MKL(替代 LAPACK 为 Armadillo 提供矩阵分解计算)。

项目使用 Visual Studio Ultimate 2012 建立,不过核心代码只有一个 cpp 文件。

全部代码托管在github.com/johnhany/QR-PCA-FaceRec


/************************************************************************/
/*                   QR-PCA-FaceRec  by John Hany                       */
/*																		*/
/*	A face recognition algorithm using QR based PCA.              		*/
/*																		*/
/*	Released under MIT license.											*/
/*																		*/
/*	Contact me at johnhany@163.com										*/
/*																		*/
/*	Welcome to my blog http://johnhany.net/, if you can read Chinese:)	*/
/*																		*/
/************************************************************************/

#include <opencv2/opencv.hpp>
#include <armadillo>
#include <iostream>
using namespace std;

int component_num = 7;

string orl_path = "G:\\Datasets\\orl_faces";

enum distance_type {ECULIDEAN = 0, MANHATTAN, MAHALANOBIS};
//double distance_criterion[3] = { 10.0, 30.0, 3.0};
double distance_criterion[3] = { 1000.0, 1000.0, 1000.0};

bool compDistance(pair<int, double> a, pair<int, double> b);
double calcuDistance(const arma::vec vec1, const arma::vec vec2, distance_type dis_type);
double calcuDistance(const arma::vec vec1, const arma::vec vec2, const arma::mat cov2, distance_type dis_type);

int main(int argc, const char *argv[]) {
	
	int class_num = 40;
	int sample_num = 10;

	int img_cols = 92;
	int img_rows = 112;
	cv::Size sample_size(img_cols, img_rows);

	arma::mat mat_sample(img_rows*img_cols, sample_num*class_num);

	//Load samples in one matrix `mat_sample`.

	for(int class_idx = 0; class_idx < class_num; class_idx++) {
		for(int sample_idx = 0; sample_idx < sample_num; sample_idx++) {

			string filename = orl_path + "\\s" + to_string(class_idx+1) + "\\" + to_string(sample_idx+1) + ".pgm";
			cv::Mat img_frame = cv::imread(filename, CV_LOAD_IMAGE_GRAYSCALE);
			cv::Mat img_sample;
			cv::resize(img_frame, img_sample, sample_size);

			for(int i = 0; i < img_rows; i++) {
				uchar* pframe = img_sample.ptr(i);
				for(int j = 0; j < img_cols; j++) {
					mat_sample(i*img_cols+j, class_idx*sample_num+sample_idx) = (double)pframe[j]/255.0;
				}
			}
		}
	}
//	cout <<	mat_sample.n_rows << endl << mat_sample.n_cols << endl << mat_sample(img_rows*img_cols/2, 0) << endl;

	//Calculate PCA transform matrix `mat_pca`.

	arma::mat H = mat_sample;
	arma::mat mean_x = arma::mean(mat_sample, 1);

	for(int j = 0; j < class_num * sample_num; j++) {
		H.col(j) -= mean_x.col(0);
	}
	H /= 1.0/sqrt(sample_num-1);

	arma::mat Q, R;
	arma::qr_econ(Q, R, H);

	arma::mat U, V;
	arma::vec d;
	arma::svd_econ(U, d, V, R.t());

//	cout << "d" << endl << d << endl;

//	arma::rowvec vec_eigen = d.head(component_num).t();
//	cout << "vec_eigen" << endl << vec_eigen << endl;

	arma::mat V_h(V.n_rows, component_num);
	if(component_num == 1) {
		V_h = V.col(0);
	}else {
		V_h = V.cols(0, component_num-1);
	}

	arma::mat mat_pca = Q * V_h;

	//Calculate eigenfaces `mat_eigen_vec`.

	arma::mat mat_eigen = mat_pca.t() * mat_sample;
//	cout << "mat_eigen" << endl << mat_eigen << endl;
	arma::mat mat_eigen_vec(component_num, class_num, arma::fill::zeros);
	vector mat_cov_list;

	for(int class_idx = 0; class_idx < class_num; class_idx++) {

		arma::vec eigen_sum(component_num, arma::fill::zeros);
		for(int sample_idx = 0; sample_idx < sample_num; sample_idx++) {
			eigen_sum += mat_eigen.col(class_idx*sample_num+sample_idx);
		}
		eigen_sum /= (double)sample_num;
		mat_eigen_vec.col(class_idx) = eigen_sum;

		mat_cov_list.push_back(arma::cov((mat_eigen.cols(class_idx*sample_num, class_idx*sample_num+sample_num-1)).t()));

//		cout << mat_cov_list[class_idx] << endl;

	}

//	cout << "mat_eigen_vec" << endl << mat_eigen_vec << endl;

/*
	cout << "dis within class" << endl;
	for(int class_idx = 0; class_idx < class_num; class_idx++) {
		for(int sample_idx = 0; sample_idx < sample_num; sample_idx++) {
			double dis = calcuDistance(mat_eigen.col(class_idx*sample_num+sample_idx), mat_eigen_vec.col(class_idx), mat_cov_list[class_idx], distance_type::MAHALANOBIS);
			cout << dis << " ";
		}
		cout << endl;
	}

	cout << "dis between classes" << endl;
	for(int class_idx = 0; class_idx < class_num; class_idx++) {
		for(int sample_idx = 0; sample_idx < class_num; sample_idx++) {
			double dis = calcuDistance(mat_eigen.col(sample_idx*sample_num), mat_eigen_vec.col(class_idx), mat_cov_list[class_idx], distance_type::MAHALANOBIS);
			cout << dis << " ";
		}
		cout << endl;
	}
*/

	//Classify new sample.

	int correct_count = 0;

	double max_dis = 0.0;

	for(int class_idx = 0; class_idx < class_num; class_idx++){
		for(int sample_idx = 0; sample_idx < sample_num; sample_idx++) {
			arma::mat mat_new_sample(img_rows*img_cols, 1);

			string filename = orl_path + "\\s" + to_string(class_idx+1) + "\\" + to_string(sample_idx+1) + ".pgm";
			cv::Mat img_new_frame = cv::imread(filename, CV_LOAD_IMAGE_GRAYSCALE);
			cv::Mat img_new_sample;
			cv::resize(img_new_frame, img_new_sample, sample_size);

			for(int i = 0; i < img_rows; i++) {
				uchar* pframe = img_new_sample.ptr(i);
				for(int j = 0; j < img_cols; j++) {
					mat_new_sample(i*img_cols+j, 0) = (double)pframe[j]/255.0;
				}
			}

			arma::mat mat_new_eigen = mat_pca.t() * mat_new_sample;

			vector<pair<int, double>> dis_list;
			for(int new_class_idx = 0; new_class_idx < class_num; new_class_idx++) {
				double dis = calcuDistance(mat_new_eigen.col(0), mat_eigen_vec.col(new_class_idx), mat_cov_list[new_class_idx], distance_type::MAHALANOBIS);
				dis_list.push_back(make_pair(new_class_idx, dis));
			}
			sort(dis_list.begin(), dis_list.end(), compDistance);

			if(dis_list[0].first == class_idx && dis_list[0].second <= distance_criterion[distance_type::MAHALANOBIS]) { 				correct_count++; 			} 			if(dis_list.back().second > max_dis) {
				max_dis = dis_list.back().second;
			}
		}
	}

	cout << "Maximum distance: " << max_dis << endl;

	double correct_ratio = (double)correct_count / (class_num * sample_num);
	cout << "Correctness ratio: " << correct_ratio * 100.0 << "%" << endl;

	cin.get();

	return 0;
}

bool compDistance(pair<int, double> a, pair<int, double> b) {
	return (a.second < b.second);
}

double calcuDistance(const arma::vec vec1, const arma::vec vec2, distance_type dis_type) {

	if(dis_type == ECULIDEAN) {
		return arma::norm(vec1-vec2, 2);
	}else if(dis_type == MANHATTAN) {
		return arma::norm(vec1-vec2, 1);
	}else if(dis_type == MAHALANOBIS) {
		arma::mat tmp = (vec1-vec2).t() * (vec1 - vec2);
		return sqrt(tmp(0,0));
	}

	return -1.0;
}

double calcuDistance(const arma::vec vec1, const arma::vec vec2, const arma::mat cov2, distance_type dis_type) {

	if(dis_type == ECULIDEAN) {
		return arma::norm(vec1-vec2, 2);
	}else if(dis_type == MANHATTAN) {
		return arma::norm(vec1-vec2, 1);
	}else if(dis_type == MAHALANOBIS) {
		arma::mat tmp = (vec1-vec2).t() * cov2.i() * (vec1 - vec2);
		return sqrt(tmp(0,0));
	}

	return -1.0;
}

分类测试

测试样本采用 The AT&T Database of Faces,原称 The ORL Database of Faces,包含取自 40 个人的样本,每人 10 幅,共 400 幅图像,图像尺寸 92*112。

at_t_datasets
样本库的链接为http://www.cl.cam.ac.uk/research/dtg/attarchive/facedatabase.html

1.欧式距离降维及分类效果:

Euclidean_Distance_result
即使将 h 设为最大的 400(样本数),其分类正确率也只能达到 99%。

2.曼哈顿距离降维及分类效果:

Manhattan_Distance_result

在 h 为 128 时,分类正确率可以达到 100%,降维能力略好于欧式距离。

3.马氏距离降维及分类效果:

Mahalanobis_Distance_result
在 h 仅仅为 7 的时候,其分类正确率就已经达到 100%了。采用马氏距离的 PCA 方法可以将人脸数据的维度从 10000 左右降到 7,降维效果非常优秀。在 h 超过 9 时,分类过程中计算的最大马氏距离超出了双精度浮点数 double 的上限。


参考文献

[1] Sharma A, Paliwal K K, Imoto S, et al. Principal component analysis using QR decomposition[J]. International Journal of Machine Learning and Cybernetics, 2013, 4(6): 679-683.

[2] Raj D. A Realtime Face Recognition system using PCA and various Distance Classifiers[J]. 2011.

[3] Turk M, Pentland A. Eigenfaces for recognition[J]. Journal of Cognitive Neuroscience, 1991, 3(1): 71-86.

 

转载自johnhany.com


Deeplearn, 版权所有丨如未注明 , 均为原创丨本网站采用BY-NC-SA协议进行授权 , 转载请注明基于 QR 和 PCA 的人脸识别
喜欢 (0)
admin
关于作者:
互联网行业码农一枚/业余铲屎官/数码影音爱好者/二次元

您必须 登录 才能发表评论!