论文下载地址：http://download.csdn.net/detail/u010411264/8720001

clear all
close all
clc
M=100;   %训练集图片数
um=100;
ustd=80;

%这一步读取与显示训练集中的图片
S=[];     %训练集矩阵
figure(1);
addpath('C:\MATLAB7\work\TrainDatabase1');  %设计raindatabase路径
for i=1:M
    str=strcat(int2str(i),'.bmp');   %连接两字符串形成图片的文件名
    eval('img=imread(str);');        %读取图片存储在img中 在matlab中执行
    subplot(ceil(sqrt(M)),ceil(sqrt(M)),i)   %把窗口分成 4 X 5 显示
    imshow(img);                     %显示读取到的图片
    if i==3
        title('训练集图像（Training）','fontsize',20)  %设计窗口标题及大少
    end
    drawnow;   %刷新
    [irow icol]=size(img);   %提取一种图片的行与列的大少
    temp=reshape(img,irow*icol,1);      %创建一个 (N1*N2)x1 的列向量
    S=[S temp];     %创建N1*N2xM 的列向量                 
end

%这部分改变所有的图像的均值和标准值,对图像进行规格化
for i=1:size(S,2)
    temp=double(S(:,i));  %获取S矩阵的第i个列矩阵并双精度化
    m=mean(temp);         %求列向量的均值
    st=std(temp);         %求temp的标准偏差 （是样本标准差）
    S(:,i)=(temp-m)*ustd/st+um;   %由此公式对S进行均值化
end

figure(2);
for i=1:M
    %str=strcat(int2str(i),'.bmp');  %读取20张图片
    img=reshape(S(:,i),irow,icol);  %把S集合中的第i列向量转为大少为irow*icol的矩阵
    %img=img';
    %eval('imwrite(img,str)');   %%执行字符串 把图片读写人当前目录
    subplot(ceil(sqrt(M)),ceil(sqrt(M)),i)  %按4X5的行列显示
    imshow(img)    %显示规格化后的图像
    drawnow;       %刷新
    if i==3
        title('规格化后的测试集图像','fontsize',20)   %窗口标题
    end
end

%这一步显示平均脸
m=mean(S,2);   %获得行的均值 (平均脸)
tmimg=uint8(m);   %converts to unsigned 8-bit integer. Values range from 0 to 255
img=reshape(tmimg,irow,icol);    %takes the N1*N2x1 vector and creates a N2xN1 matrix
figure(3);
imshow(img);   %显示平均脸
title('训练集的平均人脸','fontsize',20)  %设计平均脸窗口的标题

%这一步求出A矩阵
tmimg=reshape(img,irow*icol,1);  %把平均脸转成向量的形式
A=[];
for i=1:M
    A(:,i)=(S(:,i))-tmimg;     
end
L=A'*A;  %得到构造矩阵L
[vv dd]=eig(L);  %vv是L的特征向量构成的列向量，dd是L特征值构成的对角阵

%消去特征值为零的特征向量和特征值
v=[];
d=[];
for i=1:size(vv,2)
    if(dd(i,i)>1e-4)    %1e-4=0.0001
        v=[v vv(:,i)];
        d=[d dd(i,i)];
    end
end    %结果这循环后，吧dd变换为d  d=1xM   d为L的特征值
[B index]=sort(d);  % 对d进行从小到大排序，排序结果存在B中，index是B在d中的索引号

ind=zeros(size(index));
dtemp=zeros(size(index));
vtemp=zeros(size(v));
len=length(index);
%此循环的目的是把特征值d和特征向量v 按左右翻转了    从大到小
 for i=1:len       %循环20次
    dtemp(i)=B(len+1-i);  % 将B(len+1-i)的值给dtemp(i) 倒序了，比如 B(20)=dtemp(1) 目的是：把B（也就是L的来的特征值）临时存储在变量dtemp中  
    ind(i)=len+1-index(i);  %也是把index索引倒序了，为了和dtemp变量一样
    vtemp(:,ind(i))=v(:,i); %Vtemp存储特征向量V，按左右翻转存储
 end
 d=dtemp;  %左右翻转后把值保持好 d为L特征值按左右翻转的结果
 v=vtemp;  %左右翻转后把值保持好 v为L特征向量按左右翻转的结果

  %把矩阵L特征向量规格化
 for i=1:size(v,2)       %循环数和矩阵v的列数相等 20次
   kk=v(:,i);
   temp=sqrt(sum(kk.^2));
   v(:,i)=v(:,i)./temp;
 end
 
 %求出协方差矩阵C的特征向量
u=[];
for i=1:size(v,2)   %循环次数和矩阵v的列数相等 20次
    temp=sqrt(d(i));
    u=[u (A*v(:,i))./temp];
end 
 
% 特征向量的规格化
for i=1:size(u,2)  %循环次数为U向量个数，为20次
   kk=u(:,i);      %将U的一个向量给予kk 大少为36000X1
   temp=sqrt(sum(kk.^2)); %%取一个临时变量,命名为temp，将kk每一项元素进行平方运算.然后求和,再取平方根 赋给temp
u(:,i)=u(:,i)./temp;
end

% 显示特征脸
figure(4);
for i=1:size(u,2)  %循环次数为U向量个数，为20次
    img=reshape(u(:,i),irow,icol);  %从向量的显示转为图片大少
    img=histeq(img,255);
    subplot(ceil(sqrt(M)),ceil(sqrt(M)),i)  %分行列显示
    imshow(img)  %显示图片
    drawnow;    %刷新
    if i==3
        title('特征人脸','fontsize',20)  %设计标题
    end
end

% 找到每个人脸在训练中的权值
omega = [];   % 建立一个空矩阵omega
for h=1:size(A,2)    %循环20次
    WW=[];           %外层循环建立一个空矩阵WW
    for i=1:size(u,2)  %循环20次
        t = u(:,i)';    
        WeightOfImage = dot(t,A(:,h));
        WW = [WW; WeightOfImage];
    end
    omega = [omega WW];       %omega是M张已知人脸投影到特征子空间ui上的坐标   omega
end

%以上部分即完成了图像的读入,规格化,特征空间的训练,特征脸的形成
%并且显示出训练图像,规格化图像,均值图像和特征脸
%以下部分为识别的部分

%获得一个新图片 并对输入的图片进行规格化
InputImage = input('请输入测试图片（一张）： \n','s'); %输入一字符串 比如1.jpg
InputImage = imread(strcat('C:\MATLAB7\work\TestDatabase1\',InputImage));
figure(5)
subplot(1,2,1);      %一个窗口分1X2显示
imshow(InputImage);  %显示输入的图像
colormap('gray');    %
title('输入的测试图像','fontsize',20);  %设计窗口标题
InImage=reshape(double(InputImage),irow*icol,1);  %把输入的图片转换为(N1*N2)x1 的列向量
temp=InImage;
me=mean(temp);              %求列向量的均值，默认求列的均值
st=std(temp);               %求的标准偏差 （是样本标准差）
temp=(temp-me)*ustd/st+um;  %由此公式对图片进行均值化
NormImage = temp;           %保存规格化后的图片在NormImage变量中

Difference = temp-m;        %m是训练集合S里的均值，得到差值向量，即是(R－Ψ)
%NormImage=Difference;

%得到输入图像在ui上的坐标
p = [];                     %测试图片主成分
aa=size(u,2);               %循环20次数 是训练集图片数
for i = 1:aa
    pare = dot(u(:,i)',Difference);
    p = [p; pare];           %p存储输入图片在子空间u上的坐标                  p 
end
%重建
ReshapedImage = m + u(:,1:aa)*p;  %m为均值图像,u是特征向量,取u矩阵的1到aa列与p矩阵相乘,再与m相加
ReshapedImage = reshape(ReshapedImage,irow,icol);
%ReshapedImage = ReshapedImage';
%show the reconstructed image.
subplot(1,2,2)
imagesc(ReshapedImage);     %显示重构图片
colormap('gray');
title('通过特征向量U重建的测试图像','fontsize',20);

InImWeight=p;

ll = 1:M;
figure(6)
subplot(1,2,1)
stem(ll,InImWeight)
title('Weight of Input Face','fontsize',14)


% Find Euclidean distance
e=[];                    %测试图片与训练集的欧氏距离只差
for i=1:size(omega,2)
    q = omega(:,i);
    DiffWeight = InImWeight-q;
    mag = norm(DiffWeight);
    e = [e mag];
end

kk = 1:size(e,2);
subplot(1,2,2)
stem(kk,e)
title('Eucledian distance of input image','fontsize',14)
MaximumValue=max(e)
MinimumValue=min(e)