function [w,J] = computeLDA(data,labels)


U = unique(labels);
if length(U) ~= 2 || U(1) ~= 0 || U(2) ~= 1
    error('class labels must be 0 or 1!')
end

X1 = data(labels==1,:);
X2 = data(labels==0,:);

% number of positive and negative examples
n1 = size(X1,1);
n2 = size(X2,1);
% class means and difference of class means
me1 = mean(X1,1);
me2 = mean(X2,1);

me_diff = me1-me2;

% subtract class mean from data
%mX_plus = Xplus - repmat(m_plus, num_plus, 1);
%mX_minus = Xminus - repmat(m_minus,num_minus,1);
% compute within class scatter
S1 = (n1-1)*cov(X1);
S2 = (n2-1)*cov(X2);
tol = 1e-2;
Sw = S1+S2;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%!!!
Sw = Sw + eye(size(Sw))*mean(diag(Sw))*tol;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%!!!

%w = inv(Sw)*me_diff'; % use this for slight regularization

w = me_diff'; % use this for infinite regularization!!!
w = w/norm(w);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% test random projection:
%w = randn(size(w));
%w = w/norm(w);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Finally, compute value of the Fisher's criterion:
% J(w) = (w'*Sb*w)/(w'*Sw*w)

% Compute between class scatter matrix Sb:
Sb = me_diff'*me_diff;

% compute J(w):
J = (w'*Sb*w)/(w'*Sw*w);