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;