lda.m

K-L+LDA结合了K-L和LDA优点的一个算法

clear all;

% Load the ATT image set
%TrainX=[];
ClassClientNum=40;                                                  %类别数目
EachClassNum=10;
Class=1:40; 
P = ClassClientNum *EachClassNum;                                   %类别的号码
k = 0;
for i=1:1:ClassClientNum
    for j=1:1:EachClassNum
        image_data = imread(['C:\Program Files\MATLAB\R2006b\work\image-base\s',num2str(i),'_',num2str(j),'.bmp']);
        temp=reshape(image_data,112*92,1);        
        %TrainClass((i-1)*EachClassNum+j)=Client(i);
       % TrainX=[TrainX temp];
        k = k + 1;
        x(:,k) = image_data(:); 
        anot_name(k,:) = sprintf('%2d:%2d',i,j);   % for plot annotations
     end;
end;
nImages = k;                     %total number of images
imsize = size(image_data);       %size of image (they all should have the same size) 
nPixels = imsize(1)*imsize(2);   %number of pixels in image
x = double(x)/255;               %convert to double and normalize

avrgx = mean(x')';
for i=1:1:nImages
    x(:,i) = x(:,i) - avrgx; % substruct the average
end;
subplot(2,2,1); imshow(reshape(avrgx, imsize)); title('mean face')

%计算协方差矩阵
L = x'*x;
%计算特征值特征向量 
[eigvec,eigval] = eig(L);          
eigval=diag(eigval)';
[eigval,I]=sort(eigval);
eigval_St_t=fliplr(eigval);
eigvec_St_t=fliplr(eigvec(:,I));
% 提取前m_b个特征值和特征向量
m_b=rank(L);
eigval_St_t=eigval_St_t(:,1:m_b);
eigvec_St_t=eigvec_St_t(:,1:m_b);
%计算St的特征向量
eigvec_St=x*eigvec_St_t;

%计算子空间V
Db=diag(eigval_St_t);
V=eigvec_St*(inv(sqrt(Db)));

subplot(2,2,2); imshow(ScaleImage(reshape(V(:,1),imsize))); title('1st eigen face');
subplot(2,2,3); imshow(ScaleImage(reshape(V(:,2),imsize))); title('2st eigen face');
subplot(2,2,4); plot(diag(eigval)); title('Eigen values');

%%%%%%%%%%%%%%%%%%%%%%%% Projecting centered image vectors onto eigenspace
% Zi = V_PCA' * (Ti-m_database)
ProjectedImages_PCA =V'*x;

image_index = 12;      %index of face to be reconstructed
reconst = V*ProjectedImages_PCA;
diff = abs(reconst(:,image_index) - x(:,image_index));
strdiff_sum = sprintf('delta per pixel: %e',sum(sum(diff))/nPixels);
figure;
subplot(2,2,1); imshow((reshape(avrgx+reconst(:,image_index), imsize))); title('Reconstructed');
subplot(2,2,2); imshow((reshape(avrgx+x(:,image_index), imsize)));title('original');
subplot(2,2,3); imshow(ScaleImage(reshape(diff, imsize))); title(strdiff_sum);


%%%%%%%%%%%%%%%%%%%%%%%% Calculating the mean of each class in eigenspace
m_PCA = mean(ProjectedImages_PCA,2); % 在特征值中的总体平均值
m = zeros(m_b,EachClassNum); 
Sw = zeros(m_b,m_b); % 类内分散矩阵的初始值
Sb = zeros(m_b,m_b);  % 类间分散矩阵的初始值

for i = 1 : ClassClientNum
    m(:,i) = mean( ( ProjectedImages_PCA(:,((i-1)*EachClassNum+1):i*EachClassNum) ), 2 )';    
    
    S  = zeros(m_b,m_b);
    for j = ( (i-1)*EachClassNum+1 ) : ( i*EachClassNum )
        S = S + (ProjectedImages_PCA(:,j)-m(:,i))*(ProjectedImages_PCA(:,j)-m(:,i))';
    end
    
    Sw = Sw + S; % Within Scatter Matrix
    Sb = Sb + (m(:,i)-m_PCA) * (m(:,i)-m_PCA)'; % Between Scatter Matrix
end

%使用Fisher准则
[J_eig_vec, J_eig_val] = eig(Sb,Sw); % Cost函数 J = inv(Sw) * Sb

Jeigval=diag(J_eig_val)';
[Jeigval,I]=sort(Jeigval);
Jeigval=fliplr(Jeigval);
Jeigvec=fliplr(J_eig_vec(:,I));

%提取前边较大的特征值和特征向量V
for r=1:size(Jeigval,2)
    S(r)=sum(Jeigval(1:r))/sum(Jeigval);
end
r=min(find(S>0.99));
%Jeigval=Jeigval(:,1:i);
%V_Fisher=Jeigvec(:,1:i);
for i=1:r-1
   V_Fisher(:,i) = Jeigvec(:,i);
end

%%%%%%%%%%%%%%%%%%%%%%%% Projecting images onto Fisher linear space
% Yi = V_Fisher' * V_PCA' * (Ti - m_database) 
for i = 1 : ClassClientNum*EachClassNum
    ProjectedImages_Fisher(:,i) = V_Fisher' * ProjectedImages_PCA(:,i);
end

for i=1:1:nImages-2
    dist(i) = sqrt(dot(ProjectedImages_Fisher(1,:)-ProjectedImages_Fisher(i,:), ProjectedImages_Fisher(1,:)-ProjectedImages_Fisher(i,:))); %euclidean
end;
subplot(2,2,4); plot(dist,'.-'); title('euclidean distance from the first face');

figure;
show_faces = 1:1:nImages/2;
plot(ProjectedImages_Fisher(show_faces,nImages), ProjectedImages_Fisher(show_faces,nImages-1),'.'); title('Desomposition: Numbers indicate (Face:Expression)');
for i=show_faces
    name = anot_name(i,:);
    text(ProjectedImages_Fisher(i,nImages), ProjectedImages_Fisher(i,nImages-1),name,'FontSize',8);
end;