人脸识别 - EigenFace（PCA）

%Eigenface Code clc; clear all; !dir *.mat % eigenfeature.mat file contains the feature vectors % you may skip the delay of loading many files by using this file directly yload=input('If want to create a new eigenfeature.mat file, enter 1: '); if isempty(yload), yload=0; end if yload == 1, dirname=[1:25]; % directory names are s1 to s25 nsubject=length(dirname); filename=[1:10]; % 1.pgm to 10.pgm nx=112; ny=92; % size of the individual image %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Step 1. Preprocessing % A. load a file % B. low pass the image, and subsample the pixel at 3:1 ratio at each side % Note that you should come up with your own preprocessing steps. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % construct a low pass Gaussian filter with 5 x 5 mask, zero mean and variance=1 h=fspecial('gaussian',5,1); nxr=ceil(nx/3); nyr=ceil(ny/3); % size of reduced image veclen=nxr*nyr; % vector length after reshape reduced image to vector meanface=zeros(veclen,1); % averaged face image xmi=zeros(veclen,nsubject); % averaged face image of i-th subject data=[]; for i=1:length(dirname), for j=1:length(filename), fname=['s' int2str(dirname(i)) '/' int2str(filename(j)) '.pgm']; tmp=imread(fname); tmp1=filter2(h,tmp,'same'); % low pass anti-aliasing filter % sub-sample image at 3:1 ratio at each dimension and convert to vector vec=reshape(tmp1(1:3:112,1:3:92),veclen,1); % data = [data; vec' dirname(i)]; % raw data file, last column is subject id xmi(:,i)=xmi(:,i)+vec; % accumulate each subject's averaged face end end clear tmp tmp1 vec; xmi=xmi/length(filename); % averaged face of each subject meanface=mean(xmi')'; % averaged face of all subjects %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Step 2. Feature Extraction % here we uses eigenfaces approach. You may also use Fisher faces % or other approaches %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% vec=data(:,1:1178)'; % 1178x 250 label=data(:,1179); % 250 x 1 nsample=length(label); % # of feature vectors = 250 vec=vec-meanface*ones(1,nsample); % subtract mean value [u,s,v]=svd(vec,0); clear v; % singular vavlue decomposition % save some memory space making signular values a vector s=tril(ones(size(s)))*diag(s); s=s/max(s); % normalized to 1 save eigenfeature1.mat else load eigenfeature1.mat; end figure(1),clf for b=1:25 fname=['s' int2str(b) '/' int2str(1) '.pgm']; eval(['subplot(5,5,' int2str(b) '),']); imshow(fname) end figure(1),clf, for i=1:nsubject% nsubject=25, eval(['subplot(5,5,' int2str(i) '),']); imagesc(reshape(xmi(:,i),38,31)); end colormap('gray'); % now we need to choose # of basis (eigenfaces). We use a heuristic % criterion so that the largest ns singular values will be chosen. % the sum of singular values has a relation to the percentage energy % of the signal that are retained. % We let user to choose this number. figure,clf, stairs(s),xlabel('n'),ylabel('% energy') title('% energy of first n singular values'); ns=input('Choose eigenface dimension (< 250): '); eface=u(:,1:ns); clear u % keep eigen faces nfig=min(ns,25); figure,clf,imagesc(reshape(meanface,nxr,nyr)),title(['average face']); colormap('gray'); figure,clf colormap('gray') for i=1:nfig, eval(['subplot(5,5,' int2str(i) '),']); imagesc(reshape(eface(:,i),38,31)); end feature=vec'*eface; % nsample x ns feature vectors clear vec; figure;plot(feature); if size(feature,2) > 1 featuregraph = zeros(2000,2000); for k = 1:size(feature,1) if feature(k,1) > -3000 & feature(k,2) > -3000 & feature(k,1) < 3000 & feature(k,2) < 3000 xpos = round(feature(k,1)/3)+1000; ypos = round(feature(k,2)/3)+1000; featuregraph(xpos-4:xpos+4,ypos-4:ypos+4) = ceil(k/10); end end else featuregraph = zeros(1,2000); for k = 1:size(feature,1) xpos = round(feature(k,1)/2)+1000; featuregraph(xpos) = ceil(k/10); end figure; plot(featuregraph); end figure; imagesc(featuregraph); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Step 3. Testing % testing is done following testing protocols used in DARPA Ferret face % recognition contest. See references: % Phillips, P.J.; Hyeonjoon Moon; Rizvi, S.A.; Rauss, P.J. % "The FERET evaluation methodology for face-recognition algorithms," % IEEE Trans. Pattern Analysis and Machine Intelligence, vol. 22, no. 10, % pp. 1090- 1104, Oct. 2000. % NOTE TESTIING PROTOCOL MUST BE FOLLOWED EXACTLY. % You may choose your own distance measure. But the % We use closed-set identification test. We separate the data into test % set and query set. For each feature vector in the query set, a distance % between it and each of the testing set feature vector will be measured. % these distances will be ranked (smallest distance rank = 1). Then, we % will compute percentage of correct identification at top n ranks. % Since each subject has 10 poses, we rotate to use each post as test image % the remaining nine as query images. %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% idx=[1:nsample]; % original index of each feature vector for j=1:10, tidx{j}=idx(j:10:nsample); % test set indices 1 x 25 qidx{j}=setdiff(idx,tidx{j}); % query set indices 1 x 225 % calculate distance from query set to test set dmat{j}=dist(feature(qidx{j},:),feature(tidx{j},:)); % 225 x 25 [tmp,tmpr]=sort(dmat{j}'); % tmpr is 25 x 225 [tmp1,tmp2]=sort(tmpr); rank{j}=tmp2'; clear tmp tmp1 tmp2; % rank{j} is 225 x 25 end % testing results % compute cumulative face recognition test results % rank{1}:rank{N}: M x nsubject matrix, M: # of query samples, N: # of partitions % nsubject: # of test samples, each sample is a subject % dirname: name of the subjects % qidx{1}:qidx{N}: indices of query samples, % tidx{1}:tidx{N}: indices of test samples % label: labels of all samples including query and test samples % each row of rank{j} gives ranks of distance from q to each of the 25 test images. % 1 is closest and 25 is furtherest according to the distance measure. % We need to check the rank of the correct label. If the label of q is sk, then % we check the k-th entry of that row. For this, we first get the positions of these % entries: cnt=[]; for j=1:length(rank), pos=[label(qidx{j})*ones(1,nsubject)==ones(length(qidx{j}),1)*label(tidx{j})']; % consists of 1 at the column where label % of q is the same as label of test image. O elsewhere. 225 x 25 rk=sum([pos.*rank{j}]'); % a 1 x 225 row vector consisting of ranks of correct label %figure; %hist(rk,dirname); topk=hist(rk,dirname); % histogram on how many rank top k, 1 <= k <= 25, 1 x 25 cnt=[cnt;topk]; end figure,stairs(tril(ones(nsubject))*sum(cnt)'/sum(sum(cnt))*100) title('cumulative ranking over 10 way cross validation') xlabel('top rank n'); ylabel('% correct identification'); eigenresult4 = tril(ones(nsubject))*sum(cnt)'/sum(sum(cnt))*100;