% IPCA_JC_3.m   25th, May, 2008  made by J.Chen
% Candid Covariance-Free Incremental Principal Component Analysis. Juyang
% Weng, Yilu Zhang, Wey-Shiuan Hwang.

k = 70;
l = 0;
%%u = zeros(10304,2,k);
%%v = zeros(10304,2,k);
x = zeros(10304,40);
m = 200;                                                                   % the No.of Samples
m1 = 40;
n = 1;
for i = 1:40
    for j = 6:6
        a = imread(strcat('E:\FaceDataBase\ORL\s',num2str(i),'\',num2str(j),'.pgm'));   %a是112*92的矩陣                                                                        
        x(:,n) = (a(1:112*92))'; 
        x = double(x);
        n = n+1;                                                                  
    end
end

%中心化數據
%%%%%%%Xm = mean(x');                                                             %x轉置后每行代表一張圖片，按列求平均,返回一個行向量1*10304
for j = 1:m1         
    Xcter(:,j) = x(:,j)-Xm';
end
x=Xcter(:,1:m1);  

for n=201:240
    if mod(n,2)~=0        %奇數
        u(:,1,1) = x(:,n-200);
    else                  %偶數
        u(:,2,1) = x(:,n-200);
    end
    for i=1:k
        if i==n
            if mod(n,2)~=0
                v(:,1,i) = u(:,1,i);
            else
                v(:,2,i) = u(:,2,i);
            end
        else
            if mod(n,2)~=0
                v(:,1,i) = ((n-1-l)/n)*v(:,2,i)+((1+l)/n)*u(:,1,i)*u(:,1,i)'*(v(:,2,i)/norm(v(:,2,i)));      % Euclidean length
                u(:,1,i+1) = u(:,1,i)-u(:,1,i)'*(v(:,1,i)/norm(v(:,1,i)))*(v(:,1,i)/norm(v(:,1,i)));            
            else
                v(:,2,i) = ((n-1-l)/n)*v(:,1,i)+((1+l)/n)*u(:,2,i)*u(:,2,i)'*(v(:,1,i)/norm(v(:,1,i)));      % Euclidean length
                u(:,2,i+1) = u(:,2,i)-u(:,2,i)'*(v(:,2,i)/norm(v(:,2,i)))*(v(:,2,i)/norm(v(:,2,i))); 
            end
        end
    end
    n
end

P = [];
for i=1:k
    P(:,i) = v(:,2,i)/norm(v(:,2,i));
end
Y = P'*x;

accu = 0; 

%-------------------------------IPCA & Test--------------------------------
for i=1:40
    for j=7:10                                                             %讀入40 x 5 副測試圖像
        a=imread(strcat('E:\FaceDataBase\ORL\s',num2str(i),'\',num2str(j),'.pgm'));
        b=(a(1:10304))';
        b=double(b);
        bm = b - Xm';                                                      %將測試圖片中心化!!!
        t=P'*bm;                                                           %計算測試圖片在新特征空間中的坐標，是k*1階矩陣
        mdist = [];
        for s=1:m1
            mdist(s)=norm(t'-Y(:,s)');
        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/160                                                        % 輸出識別率

%-------------------------------IPCA & Test--------------------------------