2DPCA的原理推导与实现

文章首发于我的个人博客

1、前言

在基于PCA的人脸识别算法中，二维人脸图像要先转换为一维向量。这就导致人脸图像向量将达到高维空间，由于其协方差矩阵维度很高而且训练样本很少，因此很难评估协方差矩阵的准确性。虽然我们可以通过SVD的方法来简化计算，但是将二维图像转换为一维向量，丢失了图像行列的相关信息。而2DPCA是基于二维图像矩阵的，这种处理方法不需要事先把图像转成一维向量，图像的协方差矩阵可以通过原始图像矩阵直接构造出来。

3、2DPCA推导

2.1 首先定义一个线性变换方程

Y=AXY = AX Y=AX

上式中A为原始图像矩阵，X是一组正交基（特征子空间），Y是图像投影到特征子空间后的数据。

2.2 构造投影特征的协方差矩阵

SX=E(Y−EY)(Y−EY)TSx=E[AX−E(AX)][AX−E(AX)]TSx=E[(A−EA)X][(A−EA)X]TS_X = E(Y-EY)(Y-EY)^T \\ S_x = E[AX-E(AX)][AX-E(AX)]^T \\ S_x = E[(A-EA)X][(A-EA)X]^T SX=E(Y−EY)(Y−EY)TSx=E[AX−E(AX)][AX−E(AX)]TSx=E[(A−EA)X][(A−EA)X]T

2.3 构造图像协方差矩阵

Gt=E[(A−EA)T)(A−EA)]G_t = E[(A-EA)^T)(A-EA)] Gt=E[(A−EA)T)(A−EA)]

矩阵GtG_tGt叫做图像协方差矩阵。容易验证，GtG_tGt是一个nxn正定矩阵。我们可以使用训练样本来评估GtG_tGt，假设共有M个训练图像样本，第j个训练样本由Aj(j=1,2,...,M)A_j(j=1,2,...,M)Aj(j=1,2,...,M)由mxn矩阵构成，并且所有样本的平均图像为$\bar{A} 。那么，。那么，。那么，G_t$可简化为

Gt=1M∑j=1M(Aj−Aˉ)T(Aj−Aˉa)G_t = \frac{1}{M}\sum_{j=1}^{M}(A_j-\bar{A})^T(A_j-\bar{A}a) Gt=M1j=1∑M(Aj−Aˉ)T(Aj−Aˉa)

2.4 计算图像协方差矩阵的特征值和特征向量

得到图像协方差矩阵GtG_tGt后，计算GtG_tGt的特征值λi\lambda_iλi和特征向量viv_ivi。
对特征值进行排序，取前d个特征值对应的特征向量作为特征子空间。

w={v1,v2,...,vd}w = \{v_1,v_2,...,v_d\} w={v1,v2,...,vd}

2.5 特征提取（将图像映射到特征子空间）

Yk=wTAkY_k = w^TA_k Yk=wTAk

www为特征向量构成的矩阵（其维度为mxt），AkA_kAk为第k个图像样本。

V=[Y1,Y2,...,Yd]V = [Y_1,Y_2,...,Y_d] V=[Y1,Y2,...,Yd]

2.6 分类识别

d=∑j=1d∣∣Yk(i)−Yk(j)∣∣2d = \sum_{j=1}^{d}||Y_k^{(i)}-Y_k^{(j)}||_2 d=j=1∑d∣∣Yk(i)−Yk(j)∣∣2

2.7 图像重建

A^=wV\hat{A} = wV A^=wV

3、2DPCA图像重建结果

t=26是前t个特征值的和占所有特征值和的前99%取值。

代码实现

void PcaAchieve::PCA2D_achieve(){//数据集的划分按自己情况...std::string train_image_path = "../images/ORL_faces_train_3/s";int per_class_num = 8;   //每一个类别的样本数int class_nums = 40;    //类别数//图像样本矩阵  (40 x 112 x 92)std::vector<cv::Mat> all_train_image;//平均脸      （112 x 92）cv::Mat average_face(112, 92, CV_32FC1, cv::Scalar(0));for (int i = 1; i <= 40; i++) {std::string train_image_sub_path = train_image_path + std::to_string(i);cv::Mat per_image(112, 92, CV_32FC1, cv::Scalar(0));for (int j = 1; j <= per_class_num; j++) {std::string image_path = train_image_sub_path + "/" + std::to_string(j) + ".png";std::cout << image_path << std::endl;cv::Mat image = cv::imread(image_path, 0);image.convertTo(image, 5, 1.0/255);per_image += image;}cv::Mat average_image = per_image / 8;average_face += average_image;all_train_image.push_back(average_image);}average_face /= class_nums;cv::imshow("average_face", average_face);cv::waitKey(0);//计算协方差矩阵 (112 x 92) x (92 x 112)cv::Mat conv(112, 112, CV_32FC1, cv::Scalar(0));for (int i = 0; i < class_nums; i++) {//去均值脸cv::Mat decenter_image = all_train_image[i] - average_face;conv += decenter_image * cv::Mat(decenter_image.t());}cv::Mat eigen_values;cv::Mat eigen_vectors;cv::eigen(conv, eigen_values, eigen_vectors);float eigen_value_sum;for (int j = 0; j < eigen_values.rows; j++) {eigen_value_sum += eigen_values.at<float>(j, 0);std::cout << eigen_values.at<float>(j, 0) << std::endl;}float lamb = 0;int t = 0;for (; t < eigen_values.rows; t++) {lamb += eigen_values.at<float>(t, 0);std::cout << lamb / eigen_value_sum << std::endl;if (lamb / eigen_value_sum >= 0.99) {break;}}std::cout << "get_t:" << t << std::endl;t = 5;//(112 x t)cv::Mat L_eigen_vec(112, t, 5, cv::Scalar(0));for (int i = 0; i < t; i++) {L_eigen_vec(cv::Rect(i, 0, 1, 112)) += eigen_vectors(cv::Rect(i, 0, 1, 112));}cv::Mat temp(t, 112, 5, cv::Scalar(0));for (int i = 0; i < t; i++) {temp(cv::Rect(0, i, 112, 1)) += eigen_vectors(cv::Rect(0, i, 112, 1));}std::cout << "-------------";cv::Mat test_image = cv::imread("../images/ORL_faces_test_3/s1/10.png", 0);test_image.convertTo(test_image, CV_32FC1, 1.0/255);test_image -= average_face;//将样本映射到特征子空间 (t x 112) x (112 x 92) = (t x 92)cv::Mat Y_test = cv::Mat(L_eigen_vec.t()) * test_image;//样本重建 (112 x t) x (t x 92)cv::Mat rebuild_image = L_eigen_vec * Y_test;rebuild_image += average_face;cv::imshow("rebuild_image", rebuild_image);cv::waitKey(0);
}

4、2DPCA人脸识别结果

2DPCA识别人脸准确率

特征子空间维度	(112 x 1)	(112 x 2)	(112 x 3)	(112 x 5)	(112 x 10)	(112 x 26)	(112 x 40)
识别准确率(%)	41.25	58.75	70	73.75	83.75	90	91.25

代码

void recg_face_with_2dpca()
{/***人脸识别部分***///将训练集投影到特征子空间std::vector<cv::Mat> train_sub;for (auto iter = all_train_image.begin(); iter != all_train_image.end(); iter++) {train_sub.push_back(cv::Mat(L_eigen_vec.t()) * (*iter-average_face));}//读取测试集图像std::string test_image_path = "../images/ORL_faces_test_3/s";int right = 0;for (int i = 1; i <= 40; i++) {std::string test_image_sub_path = test_image_path + std::to_string(i);cv::Mat per_image(112, 92, CV_32FC1, cv::Scalar(0));for (int j = 0; j < 2 ; j++) {std::string image_path = test_image_sub_path + "/" + std::to_string(j+9) + ".png";cv::Mat test_image = cv::imread(image_path, 0);test_image.convertTo(test_image, 5, 1.0 / 255);test_image -= average_face;cv::Mat test_image_sub = cv::Mat(L_eigen_vec.t()) * test_image;int result = compute_face_dist(train_sub, test_image_sub);std::cout << "result:" << result << std::endl;if (result == i)right++;}}std::cout << t << "  recg right :" << right << "recg accracy:" << (float)(right / 80.0) << std::endl;
}int PcaAchieve::compute_face_dist(std::vector<cv::Mat>& train_sub, cv::Mat test_image) {std::vector<float> dist_vector;for (auto iter = train_sub.begin(); iter != train_sub.end(); iter++) {float temp = cv::norm(*iter, test_image);dist_vector.push_back(temp);}auto position = std::min_element(dist_vector.begin(), dist_vector.end());int result = position - dist_vector.begin();return result + 1;
}

5、参考资料

[1] Jian Yang, D. Zhang, A. F. Frangi and Jing-yu Yang, “Two-dimensional PCA: a new approach to appearance-based face representation and recognition,” in IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 26, no. 1, pp. 131-137, Jan. 2004, doi: 10.1109/TPAMI.2004.1261097.

[2] 图像的 2DPCA 与 2D2DPCA 特征提取

[3] 【python】双向二维PCA（2D-2D PCA）算法实现