eldadesign.m

PCA、LDA人脸检测

function LWt=ELDADesign(n,m,sampleset,classnum,samplelabel,eignum)
% ELDADESIGN will return ELDA(efficient LDA) transformation matrix
% n denotes the problem's dimension
% m is the number of samples
% sampleset stores the sample data i.e.sampleset(:,j) represents the j-th
% sample
% classnum denotes the number of classes
% samplelabel(1*m) specify the class number to which the correspond sample belongs
% eignum specify the number of discriminant vectors 

if nargout>1
    error('Too many output arguments.');
end
if nargin~=6
    error('Wrong number of input arguments.');
end
[cn,cm]=size(sampleset);
if cm~=m | cn~=n
    error('Wrong input data.');
end
[cn,cm]=size(samplelabel);
if cn~=1 | cm~=m
    error('Wrong input data');
end

%A,B分别存放每个类的平均值和所包含的样本的数目
A=zeros(n,classnum);
B=zeros(1,classnum);

for i=1:classnum
    index=find(samplelabel==i);
    for j=index
        A(:,i)=A(:,i)+sampleset(:,j);
        B(1,i)=B(1,i)+1;
    end
    A(:,i)=A(:,i)/B(1,i);
end

%求所有样本整体的平均值
average=zeros(n,1);
for j=1:m
    average=average+sampleset(:,j);
end
average=average/m;

%求解Sb矩阵
A1=A;
for i=1:classnum
    A1(:,i)=A1(:,i)-average;
end
B1=sqrt(B);
for i=1:classnum
    A1(:,i)=A1(:,i)*B1(1,i);
end
Sb=A1*A1';
Sb=Sb/m;

%求解St矩阵
for j=1:m
  sampleset(:,j)=sampleset(:,j)-average;
end
St=sampleset*sampleset';
St=St/m;

LWt1=zeros(eignum,n);

%%计算St的零子空间
St0=null(St);
St0=orth(St0);
%%计算St0的正交补空间
St0com=null(St0');
St0com=orth(St0com);

%计算第一个鉴别向量
[cm,n1]=size(St0com);
if n1==n
    [v,d]=eig(inv(St)*Sb);
    
    mark=1;
    max=d(1,1);
    for j=2:n1
        if d(j,j)>max
            max=d(j,j);
            mark=j;
        end
    end
    LWt1(1,:)=(v(:,mark))';
else
    [v,d]=eig(inv(St0com'*St*St0com)*(St0com'*Sb*St0com));
    
    mark=1;
    max=d(1,1);
    for j=2:n1
        if d(j,j)>max
            max=d(j,j);
            mark=j;
        end
    end
    LWt1(1,:)=(St0com*v(:,mark))'/sqrt(dot(St0com*v(:,mark),St0com*v(:,mark)));
end

%计算其余的鉴别向量
for i=2:eignum
    St0=[St0,(LWt1(i-1,:))'];
    St0com=null(St0');
    St0com=orth(St0com);
    
    [v,d]=eig(inv(St0com'*St*St0com)*(St0com'*Sb*St0com));
    
    mark=1;
    max=d(1,1);
    for j=2:n1-i+1
        if d(j,j)>max
            max=d(j,j);
            mark=j;
        end
    end
    LWt1(i,:)=(St0com*v(:,mark))'/sqrt(dot(St0com*v(:,mark),St0com*v(:,mark)));
end

%如果没有指定输出参数就将变换矩阵写入文件
if nargout==1
    LWt=LWt1;
else
    fid=fopen('eldawt','w');
    count=fwrite(fid,eignum,'short');
    count=fwrite(fid,LWt1,'float');
    if count~=eignum*n
        error('file write error');
    end
    fclose(fid);
end