runfacerecog_vv.m

基于PCA和LDA的人脸识别。％ Main function that reads in files, applies pre-processing, calls face ％ recognitio

save workspace_dump
disp('done')


% =======================================================================
function [P,train] = computeEigenfaces(train,Mp,plots)
% Outputs:
%   P       => eigenfaces
%   train.wt=> weights for each training face
%   train.mean => average image

% Find the image size, [Nx Ny], and the number of training images, M
[Nx Ny M] = size(train.I);

if Mp>M
    warning(sprintf(...
        ['Can''t use more principal comp than input imgaes!\n'...
            '  -> Using %d components.'],M))
    Mp = M;
end

% Use cached version of slow eigenvalue computations if found
svdFile = fullfile(train.path,'storedSvd.mat');
foundCache = exist(svdFile,'file');
if foundCache & ~train.forceNewSvd 
    % loads stored vars
    load(svdFile)
else
	% Compute EigenFaces using "training" faces
	% learn principal components from {x1, x2, ..., xn}
	% (1) find mean face, me, and
	%     differences from means, A
	X = double(reshape(train.I,[Nx*Ny M]))./256; % 1 column per face
	me = mean(X,2);
	A = X - repmat(me,[1 M]);
    clear X
	
	% (2) covariance matrix, S = A*A' (skip computing by using svd)
	% (3) partial eigenvalue decomposition S = U'*E*U
    [U,E,V] = svd(A,0); % singular val decomp much faster
	
	% (4) get sorted eigenvalues (diag of E) and eigenvectors (U)
	eigVals = diag(E);
	eigVecs = U;
    clear U V
    
    % store cache for future runs
    save(svdFile,'eigVecs','eigVals','me','A')
end
	
% (5) P' = [u1' u2' ... um'] % pick Mp principal components
P = eigVecs(:,1:Mp);        % ouput eigenfaces
lambda = eigVals(1:Mp);     % output weights

train.mean = me;

% Project each face in training set onto eigenfaces, storing weight
train.wt = P'*A;

% Reconstruct projected faces
R = P*train.wt + repmat(train.mean,[1 M]);

% Plot average face, eigenvals
if plots.intermediateOn % >> help truesize
	figure,imshow(reshape(train.mean,[Nx Ny])),title('avg face')
    if plots.savePlotsOn, saveas(gcf,'avg_face','png'), end
    figure,plot([1:length(eigVals)], eigVals,'x-'),title('\lambda strength')
    if plots.savePlotsOn, saveas(gcf,'eigval_strength','png'), end
end

% Plot eigenfaces
if plots.intermediateOn
	I = reshape(P,[Nx Ny 1 Mp]);
	for i = 1:Mp % scale for plot
        mx = max(P(:,i));
        mi = min(P(:,i));
        I(:,:,1,i) = (I(:,:,1,i)-mi)./(mx-mi);
	end 
	figure,montage(I),title('eigenfaces'); % eigenfaces
    if plots.savePlotsOn, saveas(gcf,'eigenfaces','png'), end
end
err = sum(eigVals(Mp+1:M).^2);

% Plot reconstructed images
if plots.intermediateOn
	I = reshape(R,[Nx Ny 1 M]);
	for i = 1:M % scale for plot
        mx = max(R(:,i));
        mi = min(R(:,i));
        I(:,:,1,i) = (I(:,:,1,i)-mi)./(mx-mi);
	end 
	figure,montage(I),title('reconst training images')
    if plots.savePlotsOn, saveas(gcf,'reconst_training_images','png'), end
end


% =======================================================================
function [P,train] = computeFisherfaces(train,trainClass,plots,P1)
% Outputs:
%   P       => fisherfaces
%   train.wt=> weights for each training face
%   train.mean => average image
% Summary:
%   The idea behind this approach is the maximize the ratio of
% between-class scatter to that of within-class scatter.

% Find the image size, [Nx Ny], and the number of training images, M
[Nx Ny M] = size(train.I);

% Find mean face, me, and
%     differences from means, A
X = double(reshape(train.I,[Nx*Ny M]))./256; % 1 column per face
me = mean(X,2);
A = X - repmat(me,[1 M]);

% Use cached version of slow scattermatrix computations if found
scatFile = fullfile(train.path,'storedScat.mat');
foundCache = exist(scatFile,'file');
numOfClass = length(trainClass.num);
if foundCache & ~train.forceNewScat 
    % loads stored vars
    load(scatFile)
else
	% Calculate the between-class scatter matrix, Sb
	%       and the within-class scatter matrix,  Sw
	prod = zeros(Nx*Ny);
	Sb = zeros(Nx*Ny);
	for i = 1:numOfClass
        row = trainClass.mean{i} - me;
        prod = row * row';
        Sb = Sb + prod;
	end
	
	Sw = zeros(Nx*Ny);
	for i = 1:numOfClass
        for j = (trainClass.classStartIndex(i)):(trainClass.classEndIndex(i))
            row = X(:,j) - trainClass.mean{i};
            prod = row * row';
            Sw = Sw + prod;
        end
	end
    clear prod row
    
    % store cache for future runs
    save(scatFile,'Sb','Sw')
end

% Use PCA to project into subspace
Sbb = P1.'*Sb*P1; 
Sww = P1.'*Sw*P1;
clear Sb Sw % save memory

% Current decomposition method: (from class)
% Find generalized eigenvalues & eigenvectors using eig(A,B)
[V,D] = eig(Sbb,Sww);

% Another possible method: (from class)
% 1. Note that we only care about the direction of Sw*W on m1-m2
% 2. Guess w = Sw^-1 * (m1-m2), then iterate ???

% One more possible method: (from Duda book)
% 1. Find the eigenvalues as the roots of the characteristic
%    polynomial:  
%       det(Sb - lambda(i)*Sw) = 0
% 2. Then solve for the eigenvectors w(i) directly using:
%       (Sb - lambda(i)*Sw)*w(i) = 0

% Extract eigenvalues and sort largest to smallest
Ds = diag(D);
[tmp,ndx] = sort(abs(Ds));
ndx = flipud(ndx);
% get sorted eigenvalues (diag of E) and 
% eigenvectors (project V back into full space using P1)
eigVals = Ds(ndx);
eigVecs = P1*V(:,ndx);
clear D Ds V % save a little memory

% Only keep numOfClass-1 weights

% Only keep numOfClass-1 weights, and
% Scale to make eigenvectors normalized => sum(P(:,1).^2)==1
Mp = numOfClass-1;      
lambda = eigVals(1:Mp); % output weights
P = eigVecs(:,1:Mp);    % ouput fisherfaces
P = P./repmat(sum(P.^2).^0.5,Nx*Ny,1); % normalize

train.mean = me;

% Project each face in training set onto fisherfaces, storing weight
train.wt = P.'*A;

% Reconstruct projected faces
R = P*train.wt + repmat(train.mean,[1 M]);

% Plot average face, eigenvals
if plots.intermediateOn % >> help truesize
    figure,plot([1:length(eigVals)], eigVals,'x-'),title('\lambda strength')
    if plots.savePlotsOn, saveas(gcf,'fish_eigval_strength','png'), end
end

% Plot fisherfaces
if plots.intermediateOn
	I = reshape(P,[Nx Ny 1 Mp]);
	for i = 1:Mp % scale for plot
        mx = max(P(:,i));
        mi = min(P(:,i));
        I(:,:,1,i) = (I(:,:,1,i)-mi)./(mx-mi);
	end 
	figure,montage(I),title('fisherfaces') % fisherfaces
    if plots.savePlotsOn, saveas(gcf,'fisherfaces','png'), end
end

% Plot reconstructed images
if plots.intermediateOn
	I = reshape(R,[Nx Ny 1 M]);
	for i = 1:M % scale for plot
        mx = max(R(:,i));
        mi = min(R(:,i));
        I(:,:,1,i) = (I(:,:,1,i)-mi)./(mx-mi);
	end 
	figure,montage(I),title('reconst training images')
    if plots.savePlotsOn, saveas(gcf,'fish_reconst_trainign_images','png')
    end
    
	%figure
	%imagesc(I(:,:,1,1))
    %set(gca,'Units','pixels','Position',[100 100 3*[Ny Nx]])
	%colormap gray
end


% =======================================================================
function [recog] = classifyFaces(recog,train,P,threshFace,threshClass,plots)
% Outputs:
%   recog.classNameEst  => estimated class name
%   recog.classEst      => estimated class number
%   recog.isCorrectClass=> 1 or 0 for correct/incorrect classification

% Find the image size, [Nx Ny], and the number of training images, M
[Nx Ny M] = size(train.I);
% Find the image size, [Nx Ny], and the number of recog images, M2
[Nx Ny M2] = size(recog.I);

% Init some values
Mp = length(P(1,:));
X2 = double(reshape(recog.I,[Nx*Ny M2]))./256; % 1 column per face
A2 = X2 - repmat(train.mean,[1 M2]);
% Project each face in recog set onto eigenfaces, storing weight
recog.wt = P'*A2;
% Reconstruct projected faces
R = P*recog.wt + repmat(train.mean,[1 M2]);

% Plot reconstructed images
if plots.intermediateOn
	I = reshape(R,[Nx Ny 1 M2]);
	for i = 1:M2 % scale for plot
        mx = max(R(:,i));
        mi = min(R(:,i));
        I(:,:,1,i) = (I(:,:,1,i)-mi)./(mx-mi);
	end 
	figure,montage(I),title('reconst recog images')
    if plots.savePlotsOn, saveas(gcf,'reconst_recog_images','png'), end
end

% Find euclidian distance from each recog face to each known face
recog.euDis = zeros(M,M2);
for i = 1:M2 % each recog face
    for j = 1:M % each known face class
        recog.euDis(j,i) = sqrt(sum((recog.wt(:,i) - train.wt(:,j)).^2));
    end
end

% Classifiy to Nearest-Neighbor with two thresholds:
%   threshFace => how close a euDis has to be to any face?
%   threshWho  => how close a euDis has to be to a training face to
%                 declare match
[minDis ndx] = min(recog.euDis);
	%recog.classNameTrue % truth
	%train.classNameTrue(ndx) % estimated classification
fprintf('Results of face recognition:\n')
for i = 1:M2 % each recog face
    if minDis(i) > threshFace
        recog.classNameEst{i} = 'NonFace';
    elseif minDis(i) > threshClass
        recog.classNameEst{i} = 'UnknownFace';
    else
        recog.classNameEst{i} = train.classNameTrue{ndx(i)};
        recog.classEst(i) = ndx(i);
    end
    
    recog.isCorrectClass(i) = ...
        any(recog.classEst(i) == recog.classTrue{i});
    fprintf('\trecognized %s as %s\n',...
        recog.classNameTrue{i},recog.classNameEst{i});
end

numCorrect = sum(recog.isCorrectClass);
fprintf('\t%d of %d (%d%%) faces correctly classified\n',...
    numCorrect,M2,round(numCorrect/M2*100));