% cltest.m
% function cltest(dim, numvecperclass, classdist, numClasses, numref, maxIt);
% Sample cltest(2, 80, 8, 3, 3, 800)
% 17.01.2003, Jens Ahrensfeld
% From 
% "Color Quantization of Images", Orchard, Bouman
% IEEE Trans. on Sig. Proc., vol. 39, no. 12, pp. 2677-2690, Dec. 1991.

function [cb, idx] = clusters_ob(xs,M,tse_max)
[nc,len] = size(xs);
% Root node
tree(1) = node_root(xs, 1:len);

tse = 0;
node = tree(1);
end_nodes = 1;
num_nodes = 1;
lambdas = tree(1).lambda;
map = node.q;
for m=1:M-1,
	
    if (num_nodes == 0)
        break;
    end;
    
	[v, i] = max(lambdas);
	n = end_nodes(i);
	
	node = tree(n);
    
    L_n = node.lambda;
    num_nodes = num_nodes - 1;
    
	if(node.tse < tse_max)
		break;
	end;
	
    if L_n == 0
        continue;
    end;
        
	lambdas(i) = [];
	end_nodes(i) = [];

	[node2n, node2n1] = node_split(xs, node);

	if((node2n.N * node2n1.N) == 0)
		disp('No further split');
		continue;
	end;
	
	tree(2*n) = node2n;
	tree(2*n+1) = node2n1;
	
	L2n = tree(2*n).lambda;
	L2n1 = tree(2*n+1).lambda;
	
	end_nodes = [end_nodes 2*n 2*n+1];
	lambdas = [lambdas L2n L2n1];
	
    num_nodes = num_nodes + 2;
    
	tse = tse + tree(n).tse;
	

end;
map = [tree(end_nodes).q];
lambda = [tree(end_nodes).lambda];
len = length(lambda);
[v,i] = sort(lambda);
s = i(len:-1:1);
cnt = 1;
for s=s,
	if (isempty(tree(s).C))
		break
	 end
	idx{cnt} = tree(s).C;
	cb(:,cnt) = map(:,s);
	cnt = cnt + 1;
end;
v;

% ---------------------------------------------------------
function [node_2n, node_2n1] = node_split(xs, node)
C_n = node.C;
e_n = node.e;
q_n = node.q;

i = find(e_n'*xs(:,C_n) <= e_n'*q_n);
C_2n = C_n(i);
i = find(e_n'*xs(:,C_n) > e_n'*q_n);
C_2n1 = C_n(i);

xs_2n = xs(:,C_2n);
xs_2n1 = xs(:,C_2n1);

R_n = node.R;
m_n = node.m;
N_n = node.N;

R_2n = xs_2n*xs_2n';
m_2n = sum(xs_2n,2);
N_2n = length(C_2n);

R_2n1 = R_n - R_2n;
m_2n1 = m_n - m_2n;
N_2n1 = N_n - N_2n;

node_2n = struct('R',R_2n,'m',m_2n,'N',N_2n,'C',C_2n,'q',0,'e',0,'lambda',0, 'tse', 0);
node_2n1 = struct('R',R_2n1,'m',m_2n1,'N',N_2n1,'C',C_2n1,'q',0,'e',0,'lambda',0, 'tse', 0);

node_2n = node_stat(node_2n, xs);
node_2n1 = node_stat(node_2n1, xs);

% ---------------------------------------------------------
function out = node_root(xs, C)
xs_n = xs(:,C);
R = xs_n*xs_n';
m = sum(xs_n,2);
N = length(C);

node = struct('R',R,'m',m,'N',N,'C',C,'q',0,'e',0,'lambda',0, 'tse', 0);

out = node_stat(node, xs);

% ---------------------------------------------------------
function out = node_stat(in, xs)
xs_n = xs(:,in.C);

Cov_n = in.R - 1/in.N*in.m*in.m';
[ev,es] = eig(Cov_n);
[max_v,max_i] = max(diag(es));
e_n = ev(:,max_i);
q_n = in.m/in.N;

tse_n = norm(xs_n-repmat(q_n,1,in.N))^2;
lambda_n = sum(((xs_n - repmat(q_n,1,in.N))'*e_n).^2);
in.e = e_n;
in.q = q_n;
in.tse = tse_n;
in.lambda = lambda_n;
out = in;