PCA(Principal Component Analysis,主成分分析)是一种很常用的根据变量协方差对数据进行降维、压缩的方法。它的精髓在于尽量用最少数量的维度,尽可能精确地描述数据。
PCA对数据进行降维的过程可以用下面这个动图来解释(图片摘自http://stats.stackexchange.com/a/140579/93946):
在上图中,一组位于直角坐标系的二维样本集,沿着斜线的方向有很强的相关性。所以如果我们将直角坐标系转换到斜向,也就是让横轴沿斜线方向,纵轴垂直于斜线方向。于是,在这个新的坐标系下,数据点在横轴上分布很分散,但是在纵轴方向比较集中。如果在误差允许范围内,我们完全可以将数据点在新纵轴上的坐标全部置为0,只用新横轴上的坐标来表示点的位置。这样,就完成了对数据降维的过程(即将原始直角坐标系的2个维度减小到新坐标系的1个维度)。对更高维的情况,处理过程与之类似。
PCA人脸识别
将PCA用于人脸识别的过程如下:
1.假设有400幅尺寸为100*100的图像,构成10000*400的矩阵;
2.计算均值,令
;
3.根据定义,计算协方差矩阵;
4.计算的特征值与特征向量,取前h个最大特征值所对应的特征向量,构成矩阵
;
5.矩阵可对数据降维:
,Y是h行400列的矩阵,也就是将数据从10000维降为h维。
这种做法一个明显的缺陷在于,的维度为10000×10000,直接进行奇异值分解计算量非常大。利用QR分解,作间接的奇异值分解,可以减小计算量。
利用QR分解减小计算量
基于QR分解的PCA算法步骤如下:
1.已知,其中
为d*d,H为d*n,d代表原始数据的维数,n代表样本数,d远大于n;
2.对H作QR分解,,其中Q为d*t,R为t*n,
;
3.则,对
作奇异值分解
,其中U为n*t,V为t*t,
;
4.于是,其中
;
5.由于,所以QV可将
对角化,QV为
的特征向量矩阵,
为
的特征值矩阵;
6.选取D前h个最大对角元所对应于V中的h个列,构成t*h的矩阵,则降维矩阵
;
该过程涉及对一个很大的矩阵的QR分解,和对一个较小矩阵的奇异值分解。计算量与传统方法相比较的结果如下(图片摘自[1]):
进一步,进行人脸识别的过程如下:
1.假设有c个类别,每类包含s个样本,则n=c∗s;
2.对X计算,将Y(也称特征脸)按类别计算均值,得到c个长度为h的列向量
;
3.对于未知类别的新样本x,计算,y的长度为h;
4.计算距离,取距离最小的i作为x的类标号。
距离度量d
1.欧式距离
2.曼哈顿距离
3.马氏距离
在马氏距离中,x与y分布相同,且协方差矩阵为S。加入协方差矩阵的逆矩阵的作用是,将如下图(图片部分取自http://stats.stackexchange.com/a/62147/93946)中呈椭圆分布的数据归一化到圆形分布中,再来比较距离,可以抵消不同样本集在特征空间中的分布差异。
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<uchar>(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<arma::mat> 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<uchar>(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。
样本库的链接为http://www.cl.cam.ac.uk/research/dtg/attarchive/facedatabase.html。
1.欧式距离降维及分类效果:
即使将h设为最大的400(样本数),其分类正确率也只能达到99%。
2.曼哈顿距离降维及分类效果:
在h为128时,分类正确率可以达到100%,降维能力略好于欧式距离。
3.马氏距离降维及分类效果:
在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.
arma::vec eigen_sum(component_num, arma::fill::zeros); 这行的 arma::vec是否应该写成arma::rowvec?因为矩阵分析来看,应当为h行的特征值,而vec默认是和colvec是等价的?
楼主我验证了,是我错了hhh
您好,如果我直接对QR分解完后的R矩阵进行PCA处理是不是理论上也是可行的呢
您好!71行代码H /= 1.0 / sqrt(sample_num – 1);是不是应该为class_num*sample_num?
赞同。这里是求矩阵的方差,作者在上文将读入的图像数据进行了分组,[x1, x2, …,xn],每个xn都是一个100*100行的列向量,所以在这里应该是对这些列向量求均值,观察arma::mat mat_sample(img_rows*img_cols, sample_num*class_num);这一句代码,那么列向量的总数是就是class_num*sample_num。
您好!71行代码H /= 1.0 / sqrt(sample_num – 1);是不是应该为class_num*sample_num?
你好,我想请问一下,你这篇博客的思想是直接对整张图片进行匹配是吧。这样的话当人脸角度发生变化,或者拍照的环境变化时,结果应该会受到比较大的影响,是不是去分别比对眼睛,鼻子,嘴巴,结果会更加可信呢?
你好,文中的方法其实完全是统计意义上的分类,只有样本在外观上很接近时才能保证有很高的精度。你所说的方法属于基于特征的分类,方法也有很多,效果一般也更好。比较经典的从五官位置辅助分类的方法比如香港中文大学的DeepID(2014年):http://mmlab.ie.cuhk.edu.hk/pdf/YiSun_CVPR14.pdf,他们所采用的五官定位算法:http://www.cv-foundation.org/openaccess/content_cvpr_2013/papers/Sun_Deep_Convolutional_Network_2013_CVPR_paper.pdf。
当然,人脸姿态变化比较大时,五官的变化也很大甚至发生遮挡。最近很流行的思路是用GAN(及其变种)来训练生成网络,顺便也能获得一个高精度的分类器,比如上礼拜中科院的TP-GAN:https://arxiv.org/pdf/1704.04086.pdf
楼主,有做过Android opencv的人脸识别么?
做过,用的是opencv官方文档里常见的Haar分类器。只不过目前OpenCV Android应用开发的代码更新工作还没有进行到人脸识别的章节。
第一个动图好棒
楼主,可以转载你的这篇文章吗?
当然可以:)只要标明来源,被转载也是我的荣幸呢!
你好,我用的vs2010+opencv2.4.9,编译时出现
fatal error C1083: 无法打开包括文件:“armadillo”: No such file or directory
想问一下什么原因
你的机器上是否正确配置了armadillo呢?另外还需要安装lapack或者intel mkl
我没有用Intel MKL,而是用LAPACK为Armadillo提供矩阵分解计算。编译通过了,运行时提示warning:svd():decomposition failed,error:Mat::cols(): indices out of bounds or incorrectly used .测试时arma::svd_econ(U, d, V, R.t());返回0。请问这是什么原因。
应该是Mat下标越界导致计算失败吧,能(用断点,或vs的call stack)确定是哪一行出错的吗?
我设置断点单步调试在78行时arma::svd_econ(U, d, V, R.t());svd_econ()返回值已经为0了,
运行到89行:V_h = V.cols(0, component_num-1);程序出错
后跳转到了debug.hpp中的这一段代码template<typename T1>
arma_cold
arma_noinline
static
void
arma_stop_logic_error(const T1& x)
{
#if defined(ARMA_PRINT_ERRORS)
{
get_stream_err1() << "\nerror: " << x << std::endl;
}
#endif
throw std::logic_error( std::string(x) );
}
component_num设置为多少?它应该是大于1的
我设置7,不行
错误信息:
Unhandled exception at 0x000007FEFD82A49D in test1.exe: Microsoft C++ exception: std::logic_error at memory location 0x00000000002EDC18.
这有可能是环境配置的问题吗
数据图像的尺寸是多大的?
92×112,跟你用的是同一个人脸库啊
而且也正确读入了?是不是前面哪个步骤出了问题,用断点看一下每一行输出矩阵的尺寸是否和公式符合。实在不行,把你的代码发到我的邮箱johnhany@163.com,我来看看
我单步调试发现应该是arma::svd_econ(U, d, V, R.t());失败,导致U,V被重置,然后出错。但我不知道为什么会分解 失败,代码已经发到你的邮箱,请把我看一下
我在我电脑上装上intel的MKL,然后添加到vs2013中,就可以成功运行了。我想之前的错误应该是Armadillo自带的库分解能力问题造成的。
非常感谢之前的热情帮助
不客气,能运行起来就好!
你好,我在我的vs上运行出现跟你一样的错误,我已经下载mkl并添加到vs中,还没有运行成功,求详细的解决办法?
您好 ,您这个程序是基于arm的吗?如果是,有X64下的吗
vs项目本身就是针对x64建立的,而且代码本身是可以跨平台的,只要编译器配置得当,两种平台都是可以的
能给个qq或者微信交流吗,如果可以发我邮箱,谢谢了
抱歉回复的很晚,最近我在抓紧整理一批Android Studio 2上用OpenCV的教程,很快就会放出来!
您好,我正在进行android平台上的人脸识别开发,我用JNI的方法,需要把识别后的图像(Mat类型)存入jintArray型数组,再把jintArray传回java层显示,你有什么好办法啊
这数学知识用的深啊~ 奇异值分解什么的 矩阵里面的内容吗?
奇异值分解在一般的矩阵理论书籍中都会有介绍,可以看成是对非方阵进行的特征值分解。通俗地讲,矩阵乘以自身的转置,这个乘积矩阵的特征值开根号就是原矩阵的奇异值。
请问有木有android opencv 人脸识别的教程
可以先参考一下OpenCV Android SDK里面的样例(OpenCV-android-sdk\samples\face-detection)