实验目的

1、实现PCA算法的人脸重构，即用20,40,60,80,...,160个投影来重构图像的效果。

2、实现PCA算法的人脸识别，给出不少于三个人脸数据库上 10,20,30,...,160维的人脸识别识别率，并打印出识别率变化曲线图。

3、用PCA用来进行人脸图像降维，实现3个不同数据集多个子集（或多个类）的二维和三维空间实现数据的可视化，不同类用不同颜色表示。

4、比较PCA、2DPCA、B2PCA在各个人脸数据集上的最高识别率，比较重构图像的差异。

实验步骤：

一，基于PCA算法的人脸重构与人脸识别，以及数据可视化

1，（条件：基于AR数据集的前120人，每人取前三张图像作为训练数据，后七张作为识别数据）

原图

不同维度下的重构效果（AR001-1.tif）

维度	正确率
20	62.70%
40	68.52%
60	71.33%
80	74.27%
100	74.53%
120	75.20%
140	75.41%
160	75.89%

取用不同数量特征值和正确率基于此情况下人脸识别识别率曲线

数据可视化

二维空间三维空间

2，（条件：基于orl数据集的前40人，每人取前三张图像作为训练数据，后七张作为识别数据）

原图

不同维度下重构效果（orl1_1.bmp）

维度	正确率
20	78.92%
40	83.57%
60	85.00%
80	85.71%
100	86.42%
120	86.07%
140	86.07%
160	86.07%

取用不同数量特征值和正确率基于此情况下人脸识别识别率曲线

数据可视化

二维空间三维空间

3，（条件：基于Yale数据集的前15人，每人取前三张图像作为训练数据，后七张作为识别数据）

原图

不同维度下重构效果(subject01_1.bmp)

维度	正确率
20	87.73%
40	86.77%
60	86.77%
80	86.77%
100	86.77%
120	86.77%
140	86.77%
160	86.77%

取用不同数量特征值和正确率基于此情况下人脸识别识别率曲线

数据可视化

二维空间三维空间

二，PCA,2DPCA,B2DPCA重构图像与识别率差异

条件：在ORL数据集上的重构效果（40个人，取前六张训练，后四张识别的情况下展示重构效果）

重构效果：

原图

PCA：

2DPCA:

B2DPCA

识别正确率（orl数据集，40个人，前六张训练，后四张识别）：

PCA:

最高正确率为96.3%

2DPCA:

最高正确率为97.5%

B2DPCA:

最高正确率为97.5%

尝试不同训练图像与识别图像比例后结果如下：

训练情况	训练3，识别7	训练4，识别6	训练5，识别5	训练6，识别4	训练7，识别3
PCA	86.4%	89.6%	90.0%	96.3%	96.7%
2DPCA	93.2%	93.3%	93.5%	97.5%	97.5%
B2DPCA	93.2%	93.3%	93.5%	97.5%	97.5%

从表中可以知道，2DPCA和B2DPCA在ORL训练集上的最高正确率基本相同，高于PCA的正确率，同时随着训练数据的增加，PCA与它们识别正确率的差距逐渐减小。

进一步发现，在训练图像6张，识别图像四张的情况下，PCA，2DPCA,B2DPCA识别时间如下：

PCA	2DPCA	B2DPCA
4.09s	1.59s	1.22s

可以看出，2DPCA和B2DPCA虽然在正确率上差别不大，但B2DPCA却比2DPCA更加节省时间，二者所耗时间均远远小于传统PCA算法。

代码：

PCA:

%1，测试与训练图像的录入（附加计时以查看效率）
tic
reshaped_faces=[];
%训练图像读取（以orl为例）
for i=1:40for j=1:10 a=imread(strcat('C:\Users\18023\Desktop\粤海小清华\homework\计导文件\人脸数据库\ORL56_46\orl',num2str(i),'_',num2str(j),'.bmp'));     b = reshape(a,2576,1); %将每一张人脸拉成列向量b=double(b); %uint-8转换为double类型，避免人脸展示时全灰的影响       reshaped_faces=[reshaped_faces, b];  %储存每一列转化后的列向量到矩阵reshaped_faces中end
end%用每个人前六张作为训练图，后四张作为测试图
%初始化测试集集矩阵
test_data_index = [];
%初始化训练集矩阵
train_data_index = [];  
for k=0:(size(reshaped_faces,2)/j)-1train_data_index= [train_data_index j*k+1:j*k+3];%建立一个用于查找所需训练图像的目录test_data_index= [test_data_index j*k+4:j*(k+1)];%同理，是测试图像的目录
end
train_data = reshaped_faces(:,train_data_index);%根据index读取矩阵 reshaped_faces的数据添加到矩阵train_data中
test_data = reshaped_faces(:, test_data_index);%根据index1读取矩阵 reshaped_faces的数据添加到矩阵test_data中% 2-特征向量矩阵（all_eigen_face)求得过程
% 求平均脸
mean_face = mean(train_data,2);%用mean函数求训练集矩阵的平均值,2是对每一行求平均，最后得到一个和原先尺寸一样的一个人脸，每个像素是平均颜色
% 中心化
normalizeddata= (train_data - mean_face);%得到中心化人脸（此处减去一列等于traindata每一列都减去这一列） 
%PCA特征值分解的矩阵为去中心化的矩阵乘其转置矩阵
cov= normalizeddata * normalizeddata';
%特征值分解
[eigen_vectors, eigenvalues] = eig(cov);
% 从对角矩阵获取特征值
eigen_values = diag(eigenvalues);
% 对特征值按索引进行从大到小排序
[sorted_eigen_values, index] = sort(eigen_values, 'descend'); %此处建立index获取特征值的位置，以便寻找对应的特征向量
%按照特征值的大小排列特征向量
sorted_eigen_vectors = eigen_vectors(:, index);
all_eigen_faces = sorted_eigen_vectors;%2重构过程
single_face = normalizeddata(:,1);%代码中取第一张脸进行重构
index = 1;
figure;
for dim=[20,40,60,80,100,120,140,160]eigen_faces = all_eigen_faces(:,1:dim);%降维与再次重构rebuild_face = eigen_faces * (eigen_faces' * single_face) + mean_face;subplot(2,4, index); index = index + 1;fig = imshow(mat2gray(reshape(rebuild_face, [56, 46])));title(sprintf("维度=%d", dim));    if (dim == 160)waitfor(fig);end
end% 3.人脸识别
index = 1;       
k = length(10:10:160);  %选取适当间隔
Y = zeros(1, k); 
% 遍历每个维度
for i = 1:kd = 10 * i;  % 通过i确定当前维度% 取出相应数量特征向量eigen_faces = all_eigen_faces(:, 1:d);% 测试、训练数据降维projected_train_data = eigen_faces' * (train_data - mean_face);projected_test_data = eigen_faces' * (test_data - mean_face);% 测试前准备过程correct_predict_number = 0;test_face_number = size(projected_test_data, 2);% 读取测试数据与初始化for each_test_face_index = 1:test_face_numbereach_test_face = projected_test_data(:, each_test_face_index);% 初始化最小距离min_distance = inf;predicted_label = -1;%在欧氏距离下与每一个降维后的数据进行比较for each_train_face_index = 1:size(projected_train_data, 2)% 计算距离distance = norm(each_test_face - projected_train_data(:, each_train_face_index));% 如果当前距离小于之前的最小距离，则更新最小距离和标签（1NN，取最近的一个）if distance < min_distancemin_distance = distance;predicted_label = floor((train_data_index(1, each_train_face_index) - 1) / j) + 1;endend% 获取该数据真实的标签real_label = floor((test_data_index(1, each_test_face_index) - 1) / j) + 1;% 判断预测结果是否正确if predicted_label == real_labelcorrect_predict_number = correct_predict_number + 1; % 如果预测正确，增加正确数endend % 计算正确率correct_rate = correct_predict_number / test_face_number;% 将当前维度对应的正确率存入YY(i) = correct_rate;  % 使用索引i来存储正确率，确保Y与d一致% 输出当前实验结果fprintf("维度=%d，总测试样本：%d，正确数:%d，正确率：%1f\n", d, test_face_number, correct_predict_number, correct_rate);
end
toc% 4，制作折线图使正确率随维度变化可视化
figure;
plot(10:10:160, Y, '-o');  % x轴为维度数，y轴为正确率
xlabel('维度数');
ylabel('正确率');
title('人脸识别正确率与维度变化');
waitfor(fig);%5，散点图绘制
for t=[2 3]eigen_faces=all_eigen_faces(:,1:t);projected_test_data=eigen_faces'*(test_data - mean_face);color=[];for x=1:size(test_data,2)color=[color floor((x-1)/19*2)];endif (t == 2)waitfor(scatter(projected_test_data(1, :), projected_test_data(2, :), [], color));elsewaitfor(scatter3(projected_test_data(1, :), projected_test_data(2, :), projected_test_data(3, :), [], color));endend

2DPCA

%1，测试与训练图像的录入
tic
reshaped_faces=[];
%训练图像读取（以orl为例）
for i=1:40for j=1:10%a=imread(strcat('C:\Users\18023\Desktop\粤海小清华\homework\计导文件\人脸数据库\ORL56_46\orl',num2str(i),'_',num2str(j),'.bmp'));     b = reshape(a,2576,1); b=double(b);       reshaped_faces=[reshaped_faces, b]; end
end%用每个人前六张作为训练图，后四张作为测试图
%定义测试集集矩阵
test_data_index = [];
%定义训练集矩阵
train_data_index = [];  
for k=0:(size(reshaped_faces,2)/j)-1train_data_index= [train_data_index j*k+1:j*k+6];%建立一个用于查找所需训练图像的目录test_data_index= [test_data_index j*k+7:j*(k+1)];%同理，是测试图像的目录
end
train_data = reshaped_faces(:,train_data_index);%读取矩阵 reshaped_faces的数据添加到矩阵train_data中
test_data = reshaped_faces(:, test_data_index);%读取矩阵 reshaped_faces的数据添加到矩阵test_data中
mean_face = mean(train_data,2);%2，2DPCA矩阵求取
%对于训练矩阵重新还原成二维叠加的状态
train_matrices = cell(size(train_data, 2), 1);
for i = 1:size(train_data, 2)train_matrices{i} = reshape(train_data(:, i), [56, 46]);
end
%所有训练图像矩阵的均值
mean_face_matrix = zeros(56, 46);
for i = 1:length(train_matrices)mean_face_matrix = mean_face_matrix + train_matrices{i};
end
mean_face_matrix = mean_face_matrix / length(train_matrices);% 初始化2DPCA协方差矩阵
cov_2dpca = zeros(46, 46);  
% 计算2DPCA的协方差矩阵
for i = 1:length(train_matrices)centered_matrix = train_matrices{i} - mean_face_matrix;%每一个二维数据进行去中心化cov_contribution = centered_matrix' * centered_matrix;%求取其协方差矩阵cov_2dpca = cov_2dpca + cov_contribution;%每一个二维数据的协方差矩阵，都加到最后求解的矩阵上
end
cov_2dpca = cov_2dpca / length(train_matrices);%除以所有参与的二维矩阵数，求得均值
%cov_2dpca即为最终进行协方差分解的矩阵
%特征值分解与按照特征值从大到小排列
[eigen_vectors, eigenvalues] = eig(cov_2dpca);
eigen_values = diag(eigenvalues);
[sorted_eigen_values, index] = sort(eigen_values, 'descend');
sorted_eigen_vectors = eigen_vectors(:, index);
all_eigen_faces = sorted_eigen_vectors;figure
% 3，计算重构过程，展示第一张训练图像的重构
for idx = 1:4dim = idx * 10;eigen_faces = all_eigen_faces(:, 1:dim);  % 选择前dim个特征向量% 计算投影系数test_image = train_matrices{1};  % 获取第一张训练图像centered_image = test_image - mean_face_matrix;  % 中心化图像projection = centered_image * eigen_faces;  % 投影到特征空间% 重构图像rebuild_face = mean_face_matrix + (projection * eigen_faces');  % 重构% 显示重构结果subplot(1, 4, idx);imshow(mat2gray(rebuild_face));  % 显示重构图像title(sprintf('维度=%d', dim));
end% 4.人脸识别
index = 1;       
k = length(10:10:40);
Y = zeros(1, k);
% 将测试图像转换为矩阵形式
test_matrices = cell(size(test_data, 2), 1);
for i = 1:size(test_data, 2)test_matrices{i} = reshape(test_data(:, i), [56, 46]);
end
% 遍历每个维度
for i = 1:kd = 10* i;% 取出相应数量特征向量projection_vectors = all_eigen_faces(:, 1:d);% 计算训练和测试图像的特征train_features = zeros(size(train_data, 2), 56*d);test_features = zeros(size(test_data, 2), 56*d);% 对每张训练图像进行投影for j = 1:length(train_matrices)% 计算投影projected = train_matrices{j} * projection_vectors; % 56×d% 将投影结果展平为特征向量train_features(j, :) = projected(:)';end% 对每张测试图像进行投影for j = 1:length(test_matrices)projected = test_matrices{j} * projection_vectors;test_features(j, :) = projected(:)';end% 测试前准备过程correct_predict_number = 0;test_face_number = size(test_features, 1);% 遍历每一个待测试人脸for each_test_face_index = 1:test_face_numbereach_test_face = test_features(each_test_face_index, :)'; % 读取测试用的脸特征% 初始化最小距离min_distance = inf;predicted_label = -1;% 对比每一张训练脸for each_train_face_index = 1:size(train_features, 1)% 计算距离distance = norm(each_test_face - train_features(each_train_face_index, :)');% 如果当前距离小于之前的最小距离，则更新最小距离和标签if distance < min_distancemin_distance = distance;predicted_label = floor((train_data_index(1, each_train_face_index) - 1) / j) + 1;endend% 获取实际标签real_label = floor((test_data_index(1, each_test_face_index) - 1) / j) + 1;% 判断预测结果是否正确if predicted_label == real_labelcorrect_predict_number = correct_predict_number + 1;endend% 计算正确率correct_rate = correct_predict_number / test_face_number;Y(i) = correct_rate;fprintf("维度=%d，总测试样本：%d，正确数:%d，正确率：%1f\n", d, test_face_number, correct_predict_number, correct_rate);
end
toc% 制作折线图使正确率随维度变化可视化
figure;
plot(10:10:40, Y, '-o');  % x轴为维度数，y轴为正确率
xlabel('维度数');
ylabel('正确率');
title('人脸识别正确率与维度变化');

B2DPCA

%1，测试与训练图像的录入
tic
reshaped_faces=[];
for i=1:40for j=1:10 a=imread(strcat('C:\Users\18023\Desktop\粤海小清华\homework\计导文件\人脸数据库\ORL56_46\orl',num2str(i),'_',num2str(j),'.bmp'));     b = reshape(a,2576,1); b=double(b);     reshaped_faces=[reshaped_faces, b];  end
end%用每个人前六张作为训练图，后四张作为测试图
%定义测试集集矩阵
test_data_index = [];
%定义训练集矩阵
train_data_index = [];  
for k=0:(size(reshaped_faces,2)/j)-1train_data_index= [train_data_index j*k+1:j*k+6];%建立一个用于查找所需训练图像的目录test_data_index= [test_data_index j*k+7:j*(k+1)];%同理，是测试图像的目录
end
train_data = reshaped_faces(:,train_data_index);%读取矩阵 reshaped_faces的数据添加到矩阵train_data中
test_data = reshaped_faces(:, test_data_index);%读取矩阵 reshaped_faces的数据添加到矩阵test_data中
mean_face = mean(train_data,2);%2，协方差矩阵求取
%对于训练矩阵重新还原成二维叠加的状态
train_matrices = cell(size(train_data, 2), 1);
for i = 1:size(train_data, 2)train_matrices{i} = reshape(train_data(:, i), [56, 46]);
end%求所有训练图像矩阵的均值
mean_face_matrix = zeros(56, 46);
for i = 1:length(train_matrices)mean_face_matrix = mean_face_matrix + train_matrices{i};
end
mean_face_matrix = mean_face_matrix / length(train_matrices);% 行和列方向特征值分解
%行方向
cov_2dpca_col = zeros(46, 46);
for i = 1:length(train_matrices)centered_matrix = train_matrices{i} - mean_face_matrix;%求得每一个二维数据去中心化的结果cov_contribution = centered_matrix' * centered_matrix;%行方向协方差矩阵求解，用转置乘以原矩阵cov_2dpca_col = cov_2dpca_col + cov_contribution;
end
cov_2dpca_col = cov_2dpca_col / length(train_matrices);
%列方向
cov_2dpca_row = zeros(56, 56);
for i = 1:length(train_matrices)centered_matrix = train_matrices{i} - mean_face_matrix;%求得每一个二维数据去中心化的结果cov_contribution = centered_matrix * centered_matrix';%列方向协方差矩阵求解，用原矩阵乘以其转置cov_2dpca_row = cov_2dpca_row + cov_contribution;
end
cov_2dpca_row = cov_2dpca_row / length(train_matrices);% 列方向特征向量求取
[eigen_vectors_col, eigenvalues_col] = eig(cov_2dpca_col);
eigen_values_col = diag(eigenvalues_col);
[sorted_eigen_values_col, index_col] = sort(eigen_values_col, 'descend');
sorted_eigen_vectors_col = eigen_vectors_col(:, index_col);% 行方向特征向量求取
[eigen_vectors_row, eigenvalues_row] = eig(cov_2dpca_row);
eigen_values_row = diag(eigenvalues_row);
[sorted_eigen_values_row, index_row] = sort(eigen_values_row, 'descend');
sorted_eigen_vectors_row = eigen_vectors_row(:, index_row);%3，人脸重构查看效果
figure;
for idx = 1:8dim_col = idx * 5;  % 列方向降维dim_row = idx * 5;  % 行方向降维% 行和列方向特征向量选取eigen_faces_col = sorted_eigen_vectors_col(:, 1:dim_col);eigen_faces_row = sorted_eigen_vectors_row(:, 1:dim_row);% 获取第一张训练图像并重构test_image = train_matrices{1};centered_image = test_image - mean_face_matrix;% B2DPCA投影和重构projection = eigen_faces_row' * centered_image * eigen_faces_col;rebuild_face = mean_face_matrix + eigen_faces_row * projection * eigen_faces_col';% 显示重构结果subplot(2, 4, idx);imshow(mat2gray(rebuild_face));title(sprintf('维度=%d×%d', dim_row, dim_col));
end%4,人脸识别部分
index = 1;       
k = length(10:10:40);
Y = zeros(1, k);% 将测试图像转换为矩阵形式
test_matrices = cell(size(test_data, 2), 1);
for i = 1:size(test_data, 2)test_matrices{i} = reshape(test_data(:, i), [56, 46]);
end% 遍历每个维度
for i = 1:kd = 10 * i;  % 每个方向的维度% 取出特征向量projection_vectors_col = sorted_eigen_vectors_col(:, 1:d);projection_vectors_row = sorted_eigen_vectors_row(:, 1:d);% 计算训练和测试图像的特征train_features = zeros(size(train_data, 2), d * d);test_features = zeros(size(test_data, 2), d * d);% 对每张训练图像进行双向投影for j = 1:length(train_matrices)centered_matrix = train_matrices{j} - mean_face_matrix;projected = projection_vectors_row' * centered_matrix * projection_vectors_col;train_features(j, :) = projected(:)';end% 对每张测试图像进行双向投影for j = 1:length(test_matrices)centered_matrix = test_matrices{j} - mean_face_matrix;projected = projection_vectors_row' * centered_matrix * projection_vectors_col;test_features(j, :) = projected(:)';end% 测试前准备过程correct_predict_number = 0;test_face_number = size(test_features, 1);% 遍历每一个待测试人脸for each_test_face_index = 1:test_face_numbereach_test_face = test_features(each_test_face_index, :)'; % 读取测试用的脸特征% 初始化最小距离min_distance = inf;predicted_label = -1;% 对比每一张训练脸for each_train_face_index = 1:size(train_features, 1)% 计算距离distance = norm(each_test_face - train_features(each_train_face_index, :)');% 如果当前距离小于之前的最小距离，则更新最小距离和标签if distance < min_distancemin_distance = distance;predicted_label = floor((train_data_index(1, each_train_face_index) - 1) / j) + 1;endend% 获取实际标签real_label = floor((test_data_index(1, each_test_face_index) - 1) / j) + 1;% 判断预测结果是否正确if predicted_label == real_labelcorrect_predict_number = correct_predict_number + 1;endend% 计算正确率correct_rate = correct_predict_number / test_face_number;Y(i) = correct_rate;fprintf("维度=%d，总测试样本：%d，正确数:%d，正确率：%1f\n", d, test_face_number, correct_predict_number, correct_rate);
end
toc% 制作折线图使正确率随维度变化可视化
figure;
plot(10:10:40, Y, '-o');  % x轴为维度数，y轴为正确率
xlabel('维度数');
ylabel('正确率');
title('人脸识别正确率与维度变化');