function [ind, tree_array, nodes]=split_vectors_tree(Mat, indices, level,tree_array, debug_flag, nodes, minSamples)
if length(indices) <= minSamples,
   len=length(indices);
   ind = length(tree_array)+1;
   tree_array=[tree_array len+1];
   tree_array=[tree_array indices'];
   nodes=nodes+1;
   return;
end
imgMat1=Mat(indices,2:size(Mat,2));
varVec=var(imgMat1);
%diffVec=max(imgMat1) - min(imgMat1);
%diffVec = diffVec.^2+0.01;
%varVec = varVec./diffVec;
[maxVar splitInd]=max(varVec);
x1=median(imgMat1(:,splitInd));
LR1=find(imgMat1(:,splitInd)>x1);
if length(LR1) > 0, 
   LL1=find(imgMat1(:,splitInd)<=x1);
   thresh=(min(imgMat1(LR1,splitInd))+max(imgMat1(LL1,splitInd)))/2;
else
   LL1=find(imgMat1(:,splitInd)<=x1);
   thresh=(min(imgMat1(LL1,splitInd))+max(imgMat1(LL1,splitInd)))/2;
end
LR1=find(imgMat1(:,splitInd)>thresh);
LL1=find(imgMat1(:,splitInd)<=thresh);
for i=1:length(LR1),
   LR1(i)=indices(LR1(i));
end
for i=1:length(LL1),
   LL1(i)=indices(LL1(i));
end
splitInd = splitInd + 1;
if debug_flag,
   disp([sprintf('Level %d: Splitting along dimension %d with threshold %f ... ', level, splitInd, thresh)]);
end
LL1num=Mat(LL1,splitInd);
LR1num=Mat(LR1,splitInd);
if debug_flag,
   XX=[LL1'; LL1num'];
   leftStr=[sprintf('%2d (%6.4f) ', XX)];
   XX=[LR1'; LR1num'];
   rightStr=[sprintf('%2d (%6.4f) ', XX)];
   disp([sprintf('\t Left: %s', leftStr)]);
   disp([sprintf('\tRight: %s', rightStr)]);
end
ind = length(tree_array)+1;
tree_array=[tree_array 0 0 0 splitInd thresh];
[left_ind, tree_array, nodes]=split_vectors_tree(Mat, LL1, level+1,tree_array, debug_flag, nodes, minSamples);
tree_array(ind+1)=left_ind;
[right_ind, tree_array, nodes]=split_vectors_tree(Mat, LR1, level+1,tree_array, debug_flag, nodes, minSamples);
tree_array(ind+2)=right_ind;
nodes=nodes+1;
return;
