PCA 人脸降维

1、从文件中读取图像数据（一共40个人，每人5张图片，图片大小为112*92，格式为pgm，每个人的图像单独存放在一个文件夹中）

function [imgRow,imgCol,FaceContainer,faceLabel]=ReadFaces(nFacesPerPerson, nPerson, bTest)
% 读入ORL人脸库的指定数目的人脸前前五张(训练)
%
% 输入：nFacesPerPerson --- 每个人需要读入的样本数，默认值为 5
%       nPerson --- 需要读入的人数，默认为全部 40 个人
%       bTest --- bool型的参数。默认为0，表示读入训练样本（前5张）；如果为1，表示读入测试样本（后5张）
%
% 输出：FaceContainer --- 向量化人脸容器，nPerson * 10304 的 2 维矩阵，每行对应一个人脸向量

if nargin==0 %default value
    nFacesPerPerson=5;%前5张用于训练
    nPerson=40;%要读入的人数（每人共10张，前5张用于训练）
    bTest = 0;
elseif nargin < 3
    bTest = 0;
end

img=imread('Data/ORL/S1/1.pgm');%为计算尺寸先读入一张
[imgRow,imgCol]=size(img);


FaceContainer = zeros(nFacesPerPerson*nPerson, imgRow*imgCol);
faceLabel = zeros(nFacesPerPerson*nPerson, 1);

% 读入训练数据
for i=1:nPerson
    i1=mod(i,10); % 个位
    i0=char(i/10);
    strPath='Data/ORL/S';
    if( i0~=0 )
        strPath=strcat(strPath,'0'+i0);
    end
    strPath=strcat(strPath,'0'+i1);
    strPath=strcat(strPath,'/');
    tempStrPath=strPath;
    for j=1:nFacesPerPerson
        strPath=tempStrPath;
        
        if bTest == 0 % 读入训练数据
            strPath = strcat(strPath, '0'+j);
        else
            strPath = strcat(strPath, num2str(5+j));
        end
        
        strPath=strcat(strPath,'.pgm');
        img=imread(strPath);
       
        %把读入的图像按列存储为行向量放入向量化人脸容器faceContainer的对应行中
        FaceContainer((i-1)*nFacesPerPerson+j, :) = img(:)';
        faceLabel((i-1)*nFacesPerPerson+j) = i;
    end % j
end % i

% 保存人脸样本矩阵
save('FaceMat.mat', 'FaceContainer')

　　2.主函数中调用上面的ReadFaces函数和fastPCA函数

function main(k)
% ORL 人脸数据集的主成分分析
%
% 输入：k --- 降至 k 维

% 定义图像高、宽的全局变量 imgRow 和 imgCol，它们在 ReadFaces 中被赋值
global imgRow;
global imgCol;

% 读入每个人的前5副图像
nPerson=40;
nFacesPerPerson = 5;
display('读入人脸数据...');
[imgRow,imgCol,FaceContainer,faceLabel]=ReadFaces(nFacesPerPerson,nPerson);
display('..............................');


nFaces=size(FaceContainer,1);%样本（人脸）数目
display('PCA降维...');
% LowDimFaces是200*20的矩阵, 每一行代表一张主成分脸(共40人，每人5张)，每个脸20个维特征
% W是分离变换矩阵, 10304*20 的矩阵
[LowDimFaces, W] = fastPCA(FaceContainer, 20); % 主成分分析PCA
visualize_pc(W);%显示主成分脸，详见第4步的图像显示
save('LowDimFaces.mat', 'LowDimFaces');
display('计算结束。');

　　3、fastPCA算法

function [pcaA V] = fastPCA( A, k )
% 快速PCA
%
% 输入：A --- 样本矩阵，每行为一个样本
%      k --- 降维至 k 维
%
% 输出：pcaA --- 降维后的 k 维样本特征向量组成的矩阵，每行一个样本，列数 k 为降维后的样本特征维数
%      V --- 主成分向量

[r c] = size(A);

% 样本均值
meanVec = mean(A);

% 计算协方差矩阵的转置 covMatT
Z = (A-repmat(meanVec, r, 1));
covMatT = Z * Z';

% 计算 covMatT 的前 k 个本征值和本征向量
[V D] = eigs(covMatT, k);

% 得到协方差矩阵 (covMatT)' 的本征向量
V = Z' * V;

% 本征向量归一化为单位本征向量
for i=1:k
    V(:,i)=V(:,i)/norm(V(:,i));
end

% 线性变换（投影）降维至 k 维
pcaA = Z * V;

% 保存变换矩阵 V 和变换原点 meanVec
save('PCA.mat', 'V', 'meanVec');

4、fastPCA函数的输出主分量阵W，是一个10340*20的矩阵，每列是一个10340维的主分量（样本协方差矩阵的本征向量），用112*92的分辨率来显示，

visualize_pc（）代码如下

function visualize_pc(E)
% 显示主成分分量（主成分脸，即变换空间中的基向量）
%
% 输入：E --- 矩阵，每一列是一个主成分分量


[size1 size2] = size(E);
global imgRow;
global imgCol;
row = imgRow;
col = imgCol;

if size2 ~= 20
   error('只用于显示 20 个主成分');
end;

figure
img = zeros(row, col);
for ii = 1:20
    img(:) = E(:, ii);
    subplot(4, 5, ii);
    imshow(img, []);
end

5、补充：基于主分量的人脸重建

　　

function [ xApprox ] = approx( x, k )
% 用 k 个主成分分量来近似（重建）样本 x
%
% 输入：x --- 原特征空间中的样本，被近似的对象
%       k --- 近似（重建）使用的主分量数目
%
% 输出：xApprox --- 样本的近似（重建）

% 读入 PCA 变换矩阵 V 和 平均脸 meanVec
load Mat/PCA.mat

nLen = length(x);

xApprox = meanVec;

for ii = 1:k
    xApprox=xApprox+((x-meanVec)*V(:,ii))*V(:,ii)';
end

function displayImage(v,row,col)
% 向量 v 为一幅图像按列存储得到的行向量，函数 displayImage(...) v 转化成原始的 row * col 图像，并显示之
%
% 输入：v --- 一幅图像按列存储得到的行向量
%      row --- 原始图像的行数
%      col --- 原始图像的列数

[m n] = size(v);


if m ~= 1
   error('v必须为行向量');
end;
if  row*col ~= n
   error('图像尺寸与输入向量 v 的维数不符');
end;

% 原始图像I
I=zeros(row,col);
I(:)=v(:);


figure;
imagesc(I);
colormap(gray);
axis image;