kfda.m

It is for Face Recognition

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function kfda()
%   z.li, 04-18-2004
%   fisher discriminant analysis with kernel mapping
% function dependency:
%   - ker(type, x1, x2)
% input:
%   K - number of  classes in data
%   x - data: d x n
%   y - label: 1 x n
%   kerType - kernel type
% output:
%   A  - the transform: K x d
%   ev  - eigen values
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%function [A,ev]=kfda(K, x, y, kerType)
function [A,ev]=kfda(K, x, y, kerType)
 dbg = 'y';
 if dbg == 'y'
   kerType = 'grbf';
   K = 2; m=32;
   styl = ['+', '.', 'x'];
   axis([0 100 0 100]); hold on;
   x = zeros(2, K*m);
   y = zeros(1, K*m);
   for j=1:K
       fprintf('\n new class: ');
       for t=1:m
          [px py]=ginput(1);
          plot(px, py, styl(j));
          fprintf('%c', styl(j));
          x(:,(j-1)*m+t) = [px py]';
          y((j-1)*m+t) = j;
       end
   end
 end

 % collect labeling
 [d, N] = size(x);          % sample dimension and total number

 % kernel sample mean for all
 M = zeros(N, 1);
 for j=1:N      % N dimension mean
   for k=1:N    % for all samples
      M(j) = M(j) + ker(kerType, x(:,j), x(:,k));
   end
 end
 M = M/N;

 mj = zeros(N, K);          % class means, N-dimensional
 nj = zeros(1, K);          % number of samples per class

 % compute Sw - intra class scatter: dxd
 Sw = zeros(N,N);
 % process all classes
 for j=1:K              % for all class j
   % get class j's samples
   indx =[];indx = find(y==j);
   nj(j) =length(indx);
   xj=zeros(d, nj(j)); xj = x(:,indx);

   % compute kernel sample mean for class j
   for k=1:N            % for all dimension k
      for t =1:nj(j)    % for all class samples
         mj(k,j) = mj(k,j)+ker(kerType, x(:,k), xj(:,t));
      end
   end
   mj(:,j) = mj(:,j)/nj(j);

   % compute kernel intra scatter
   ks = zeros(N,1);     % kernel sample
   for k=1:N            % for all dimension k
      for t=1:nj(j)     % for all class samples
         % compute kernel sample t
         for h=1:N
            ks(h) = ker(kerType, x(:,h), xj(:,t));
         end
         % intra scatter:
         Sw = Sw + (ks-mj(:,j))*(ks-mj(:,j))';
      end
   end
 end

 % compute Sb - inter class scatter: NxN
 Sb = zeros(N, N);
 for j=1:K
   Sb = Sb + nj(j)*(mj(j)-M)*(mj(j)-M)';
 end

 % find fisher discriminant subspace
 [v,d]=eig(Sb, Sw);
  
 A = v(:, 1:K-1); ev = diag(d);

 % dbg plot
 if dbg == 'y'
   hold on;
   for j=1:K
      plot(mj(1,j), mj(2,j), 'or');
   end
   plot(M(1), M(2), 'or'); plot(M(1), M(2), '+r');
 end