computeeigenfaces.m

face recognition with lda+pca

%function computeEigenfaces()
clear all;

%req='Yes';

prompt = {'Enter number of face in database :','Enter your address of classes:','Enter your number of classes:','Enter your number face in class'};
dlg_title = 'Input for peaks function';
num_lines = 1;
def = {'0','G:\projects\face detection\yalefaces\1 (','0','0'};
answer = inputdlg(prompt,dlg_title,num_lines,def);
path=cell2mat(answer(2));
count_of_face=cell2mat(answer(1));
count_of_face=str2num(count_of_face);

numberofclass=cell2mat(answer(3));
numberofclass=str2num(numberofclass);

num_face_in_class=cell2mat(answer(4));
num_face_in_class=str2num(num_face_in_class);
%path='G:\projects\face detection\yalefaces\1 (';
f=').jpg';
for i=1:count_of_face%130
    tmp=int2str(i);
    p=[path,tmp,f];
    database(:,:,:,i)=imread(p);
    [nx mx rgb ss]=size(database);
    if (rgb>2)
       database_gray(:,:,i)=rgb2gray(database(:,:,:,i));
    else 
        database_gray(:,:,i)=database(:,:,:,i);
    end
   % l(:,:,:,i)
end 


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

	%     differences from means, A
	X = double(reshape(database_gray,[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
  % (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
plots_intermediateOn=1;
plots_savePlotsOn=1;
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)

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

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

numOfClass =numberofclass; %length(trainClass.num);

 j=1;
 k=1;
 init=1;
 end_face=num_face_in_class;%2
 sum_person=0;
numberofclass_org=numberofclass;
numberofclass=count_of_face/numberofclass;

     
  for j=1:numberofclass
      
     for k=init:end_face
         if (k<55)
            sum_person =sum_person+(1/num_face_in_class)*(A(:,k));
         end
     end
       avg_person(:,j)=sum_person;
       init=init+num_face_in_class;
       end_face=end_face+num_face_in_class;%+1;
       sum_person=0;
  end

%subtract them from training face 
    j=1;
   for i=1:count_of_face
       sub(:,i)=A(:,i)-avg_person(:,j);
        if (mod(i,numberofclass_org)==0) && (i~=count_of_face) 
            j=j+1;
        end
   end
    
     k=1;
     sum_temp=0;
     init=1;
     counter=num_face_in_class;

     
    for i=1:numberofclass%13
        for k=init:counter
            if k<55
                 sum_temp=sum_temp+sub(:,k)*sub(:,k)';
            end
        end
           scat(:,:,i)=sum_temp;
           init=init+num_face_in_class;
           counter=counter+num_face_in_class;%+1;
           sum_temp=0;          
    end

%within class
    SW=0;
    for i=1:numberofclass%13
       SW=SW+scat(:,:,i);
    end
%comput between class

    SB=0;
 
    for f=1:numberofclass
         SB=SB+num_face_in_class*(avg_person(:,f)-train_mean)*(avg_person(:,f)-train_mean)';
     end
    

    
    Sb=SB;
    Sw=SW;
    
    % store cache for future runs
   % save(scatFile,'Sb','Sw')
%end
P1=P;
% Use PCA to project into subspace
ST=Sb+Sw;
Sbb = P1.'*Sb*P1; 
Sww = P1.'*Sw*P1;
STT=P1.'*ST*P1;
clear Sb Sw % save memory

% Find generalized eigenvalues & eigenvectors using eig(A,B)
[V,D] = eig(STT,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, 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
    
end

close all
% =======================================================================
prompt = {'Enter your address save train_wt :','Enter your adress save P:','Enter your address train_mean:','Enter your address face for detection'};
dlg_title = 'Input for peaks function';
num_lines = 1;
def = {'d:\hashemi\','d:\hashemi\','d:\hashemi\','d:\hashemi\1 (50).jpg'};
answer = inputdlg(prompt,dlg_title,num_lines,def);
address_train_wt=cell2mat(answer(1));
address_P=cell2mat(answer(2));
address_train_mean=cell2mat(answer(3));
path_pattern=cell2mat(answer(4));
im=imread(path_pattern);
threshFace=count_of_face;%130
threshClass=numberofclass_org;%10

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

save (address_train_wt,'train_wt');
save (address_P,'P');
save (address_train_mean,'train_mean');
%load train_mean;
%load train_wt;


im=imresize(im,[Nx Ny]);
m1=1;
[nx ny m1]=size(im);
if m1>2 
    recog_I=rgb2gray(im);%database_gray(:,:,67);%database(:,:,78);%feaD(:,:,150);%database_gray(:,:,50);%feaD(:,:,150);%un(:,:,20);
else 
    recog_I=im;
end
    [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
close all;
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

[minDis ndx] = min(recog_euDis);
	
fprintf('Results of face recognition:\n')
for i = 1:M2 % each recog face
    if minDis(i) > threshFace
        errordlg('NoN face');
    elseif minDis(i) > threshClass
         errordlg('UnknownFace');
       
    else if round(minDis)<1
        str1='Success face index =';
        str2=num2str(ndx);
        close all;
       subplot(1,2,1), subimage(imresize(recog_I,[200 200]));title('Pattern face','FontSize',8,'Color','r')
       subplot(1,2,2), subimage(imresize(database_gray(:,:,ndx),[200 200]));title('Pattern Recognition of facce','FontSize',8,'Color','g')
        message=[str1,str2];
        warndlg(message,'success');
       
        else  
           errordlg('other face');
    end
    end
end