function [I,Pkx]=bayescln(X,MI,SIGMA,Pk)
% BAYESCLN Bayes classifier for Gaussian distributiuon.
%  [I,Pkx]=bayescln(X,MI,SIGMA,Pk)
%
% This function classifies into the class according to the 
% a posteriori probability calculated using Bayes' formula 
% under the condition that the probability densities p(x|k) 
% have normal (Gaussian) distribution. Each class is described
% by one Gaussian (couple mi, sigma) p(x|k) and its a prior
% probability. 
%
% Input:
%  (notation N - dimension of the feature space, K - number of classes,
%            L - number of input vectors)
%   X [NxL]  Matrix of L input vectors (observations) to be
%              classified, X=[x_1,x2,...,x_L].
%   MI [NxK] Matrix of K vectors of mean values of p(x|k).
%     MI=[mi_1,mi_2,...,mi_K], where mi_j is a column vector [Nx1] for class j.
%   SIGMA [(NxN)xK] Matrix of covariance matrices of the density p(x|k).
%     SIGMA=[sigma_1,sigma_2,...,sigma_K], sigma_j is the covariance
%              matrix corresponding to the class j.
%   Pk [1xK] Vector with a priori probability densities. Pk(j) is an
%     a priori probability of class j.
%
% Output:
%   I [1xL] Vector of class identifiers into which the corresponding
%      points (observations) in the feature space were classified. 
%      I(j) is an identifier for vector X(:,j) and can assume integer
%      values from 1 to K.
%   Pkx [KxL] a posterior probabilities.
%
% See also LINCLASS.
%

% Statistical Pattern Recognition Toolbox, Vojtech Franc, Vaclav Hlavac
% (c) Czech Technical University Prague, http://cmp.felk.cvut.cz
% Written Vojtech Franc (diploma thesis) 19.03.2000 as bayescla
% Modifications
% 21-may-2001, V.Franc, Pkx added to putput arguments
% 25. 6.00 V. Hlavac, rewritten, renamed to bayescln,
%          comments in English.

K=size(MI,2);  % number of classes
DIM=size(X,1); % dimension
N=size(X,2);   % number of observations

%%%%%%%%%%

Pkx=zeros(N,K);
for k=1:K,
   logPk=log(Pk(k));
   halfLogDetSigma=0.5*log(det(SIGMA(:,(k-1)*DIM+1:k*DIM)));

   Pkx(:,k)=repmat(logPk,N,1)-repmat(halfLogDetSigma,N,1)-0.5*...
       mahalan(X,MI(:,k), SIGMA(:,(k-1)*DIM+1:k*DIM))';
end

Pkx=Pkx';

[tmp,I]=max(Pkx);


return;