/************************************************************************/
/* QR-PCA-FaceRec by
John Hany */
/* */
/* A face recognition algorithm using
QR based PCA. */
/* */
/* Released under MIT license. */
/* */
/* */
/* Welcome to my blog http://johnhany.net/,
if you can read Chinese:) */
/* */
/************************************************************************/
#include <opencv2/opencv.hpp>
#include <armadillo>
#include <iostream>
usingnamespacestd;
intcomponent_num=7;
stringorl_path="G:\\Datasets\\orl_faces";
enumdistance_type{ECULIDEAN=0,MANHATTAN,MAHALANOBIS};
//double distance_criterion[3] = { 10.0,
30.0, 3.0};
doubledistance_criterion[3]={1000.0,1000.0,1000.0};
boolcompDistance(pair<int,double>a,pair<int,double>b);
doublecalcuDistance(constarma::vecvec1,constarma::vecvec2,distance_typedis_type);
doublecalcuDistance(constarma::vecvec1,constarma::vecvec2,constarma::matcov2,distance_typedis_type);
intmain(intargc,constchar*argv[]){
intclass_num=40;
intsample_num=10;
intimg_cols=92;
intimg_rows=112;
cv::Sizesample_size(img_cols,img_rows);
arma::matmat_sample(img_rows*img_cols,sample_num*class_num);
//Load
samples in one matrix `mat_sample`.
for(intclass_idx=0;class_idx<class_num;class_idx++){
for(intsample_idx=0;sample_idx<sample_num;sample_idx++){
stringfilename=orl_path+"\\s"+to_string(class_idx+1)+"\\"+to_string(sample_idx+1)+".pgm";
cv::Matimg_frame=cv::imread(filename,CV_LOAD_IMAGE_GRAYSCALE);
cv::Matimg_sample;
cv::resize(img_frame,img_sample,sample_size);
for(inti=0;i<img_rows;i++){
uchar*pframe=img_sample.ptr<uchar>(i);
for(intj=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::matH=mat_sample;
arma::matmean_x=arma::mean(mat_sample,1);
for(intj=0;j<class_num
*sample_num;j++){
H.col(j)-=mean_x.col(0);
}
H/=1.0/sqrt(sample_num-1);
arma::matQ,R;
arma::qr_econ(Q,R,H);
arma::matU,V;
arma::vecd;
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::matV_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::matmat_pca=Q
*V_h;
//Calculate
eigenfaces `mat_eigen_vec`.
arma::matmat_eigen=mat_pca.t()*mat_sample;
// cout << "mat_eigen" << endl << mat_eigen
<< endl;
arma::matmat_eigen_vec(component_num,class_num,arma::fill::zeros);
vector<arma::mat>mat_cov_list;
for(intclass_idx=0;class_idx<class_num;class_idx++){
arma::veceigen_sum(component_num,arma::fill::zeros);
for(intsample_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.
intcorrect_count=0;
doublemax_dis=0.0;
for(intclass_idx=0;class_idx<class_num;class_idx++){
for(intsample_idx=0;sample_idx<sample_num;sample_idx++){
arma::matmat_new_sample(img_rows*img_cols,1);
stringfilename=orl_path+"\\s"+to_string(class_idx+1)+"\\"+to_string(sample_idx+1)+".pgm";
cv::Matimg_new_frame=cv::imread(filename,CV_LOAD_IMAGE_GRAYSCALE);
cv::Matimg_new_sample;
cv::resize(img_new_frame,img_new_sample,sample_size);
for(inti=0;i<img_rows;i++){
uchar*pframe=img_new_sample.ptr<uchar>(i);
for(intj=0;j<img_cols;j++){
mat_new_sample(i*img_cols+j,0)=(double)pframe[j]/255.0;
}
}
arma::matmat_new_eigen=mat_pca.t()*mat_new_sample;
vector<pair<int,double>>dis_list;
for(intnew_class_idx=0;new_class_idx<class_num;new_class_idx++){
doubledis=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;
doublecorrect_ratio=(double)correct_count/(class_num
*sample_num);
cout<<"Correctness
ratio: "<<correct_ratio
*100.0<<"%"<<endl;
cin.get();
return0;
}
boolcompDistance(pair<int,double>a,pair<int,double>b){
return(a.second<b.second);
}
doublecalcuDistance(constarma::vecvec1,constarma::vecvec2,distance_typedis_type){
if(dis_type==ECULIDEAN){
returnarma::norm(vec1-vec2,2);
}elseif(dis_type==MANHATTAN){
returnarma::norm(vec1-vec2,1);
}elseif(dis_type==MAHALANOBIS){
arma::mattmp=(vec1-vec2).t()*(vec1-vec2);
returnsqrt(tmp(0,0));
}
return-1.0;
}
doublecalcuDistance(constarma::vecvec1,constarma::vecvec2,constarma::matcov2,distance_typedis_type){
if(dis_type==ECULIDEAN){
returnarma::norm(vec1-vec2,2);
}elseif(dis_type==MANHATTAN){
returnarma::norm(vec1-vec2,1);
}elseif(dis_type==MAHALANOBIS){
arma::mattmp=(vec1-vec2).t()*cov2.i()*(vec1-vec2);
returnsqrt(tmp(0,0));
}
return-1.0;
}