function [th, U, y] = training(dog,cat,feature)

nd = size(dog,2);
nc = size(cat,2);

X = [dog, cat];

% SVD decomposition
[U,S,V] = svd(X,0);

U = U(:,1:feature);
animals = S*V';
dog0 = animals(1:feature, 1:nd);
cat0 = animals(1:feature, nd+1:nd+nc);

% LDA 
md = mean(dog0,2);
mc = mean(cat0,2);

Sw = 0;
for i=1:nd
    Sw = Sw + (dog0(:,i)-md)*(dog0(:,i)-md)';
end
for i=1:nc
    Sw = Sw + (cat0(:,i)-mc)*(cat0(:,i)-mc)';
end

Sb = (md-mc)*(md-mc)';

[V2,D] = eig(Sb,Sw);
[lambda, ind] = max(abs(diag(D)));
y = V2(:,ind);
y = y/norm(y);

vdog = y'*dog0;
vcat = y'*cat0;

if mean(vdog) > mean(vcat)  % dog < threshold < cat
    y = -y;
    vdog = -vdog;
    vcat = -vcat;
end

sortdog = sort(vdog);  % inceasing order
sortcat = sort(vcat);

t1 = length(sortdog);
t2 = 1;

while sortdog(t1) > sortcat(t2)
    t1 = t1 - 1;
    t2 = t2 + 1;
end

th = 0.5*(sortdog(t1)+sortcat(t2));



% showing the distribution
figure(1)
hist(sortdog,30)
hold on
plot([th th],[0,10],'r');

figure(2)
hist(sortcat,30)
hold on
plot([th th],[0,10],'r');








    
    
    

