%--------------------------------------------------------------
% 
% 
%--------------------------------------------------------------

%------------ function name ---------------- 

%Description	: 基于2DPCA的人臉識別
%Author	:	gongzzzz
%time	:	07-5-23
%Input  : 
%
%Output :
%參考文獻：	Yang Jian,Zhang David,F.Frangi Alejandro,Yang Jingyu.
%		Two-dimensional PCA:A new approach to appearance-based
%		face representation and recognition.IEEE Trans.on Pattern %
%		Analysis and Machine Intelligence,26(2004)1,131-137.
%
function [] = TwoDPCA()
	
	imgHei = 112;
	imgWid = 92;
	M = 200;	% training images number
	% 1 Get Average face from 'model.mat',if this file is 
	%   not exsit, please compute an average face first.
	load e:\ORL\model.mat;
	clear base dsort vsort allsamples xmean
	Avg = reshape(samplemean, imgHei,imgWid);

	% 2 Get 'Gt' matrix, Gt = (Ai - Avg)'((Ai - Avg))/M
	%   Get all test image
	Gt=zeros(imgWid, imgWid);
	for i=1:40
	    for j=1:5
	      a=imread(strcat('e:\ORL\s',num2str(i),'\',num2str(j),'.jpg'));
	      a=double(a);
	      Gt= Gt + (a - Avg)'*(a - Avg);
	      allsamples (:,:, (i-1) * 5 + j) = a;
	  end
	end

	Gt = Gt / M;


	% 3 Eigenvectors and Eigenvalues of Gt
	[v d]=eig(Gt);
	d1=diag(d);

	% 4 Sort descending with respect to eigenvalues
	[d2 index]=sort(d1); % descending 
	cols=size(v,2);	% 
	for i=1:cols
		vsort(:,i) = v(:, index(cols-i+1) ); % vsort 
		dsort(i)   = d1( index(cols-i+1) );  % dsort 
	end

	% 5 Recognition
	accu = 0;

	%5.1 Training image coordinates
	for (k=1: M)
		testFea(:,:,k) = allsamples(:,:,k) * vsort;
	end
	
	%5.2 Test
	for i=1:40
	    for j=6:10 %Read test images 40 x 5
		a=imread(strcat('e:\ORL\s',num2str(i),'\',num2str(j),'.jpg'));
		b=double(a);
		tcoor= b * vsort; % feature matrix of current test image
		for k=1:200 
			mdist(k)=Dis(tcoor, testFea(:,:,k));
		end;
		%三階近鄰	
		[dist,index2]=sort(mdist);
		class1=floor( (index2(1)-1)/5 )+1;
		class2=floor((index2(2)-1)/5)+1;
		class3=floor((index2(3)-1)/5)+1;
		if class1~=class2 && class2~=class3
		    class=class1;
		elseif class1==class2
		    class=class1;
		elseif class2==class3
		    class=class2;
		end;
		if class==i
		    accu=accu+1;
		end;
	    end;
	end;
	accuracy=accu/200 %輸出識別率

	% 將模型保存
	save('e:\ORL\model2D.mat', 'testFea', 'dsort', 'vsort', 'allsamples', 'Avg');	

% Distance between 'a' and 'b'
function [d] = Dis(a,b)
	
	c = a - b;
	[ir ic] = size(c);
	
	d = 0;
	for (k=1:ic)
		d = d + norm(c(:,k));
	end