function [t_best, var_res]=find_threshold(h)

var_res = -10.;
P = h/sum(sum(h));
q1 = 0;
mu1 = 0;
K = length(h);
mu = sum(P.*(1:K));
var_total = sum(P.*((1:K)-mu).^2);
for t=1:(K-1),
   mu1 = (q1*mu1+t*P(t))/(q1+P(t));
   q1 = q1+P(t);
   mu2 = (mu-q1*mu1)/(1-q1);
   cross_var = q1*(1-q1)*(mu1-mu2).^2;
   %disp([sprintf('t=%d mu1=%10.8f mu2=%10.8f q1=%10.8f cross=%10.8f',t,mu1,mu2,q1,cross_var)]);
   if var_res < cross_var,
     var_res = cross_var;
     t_best = t+1;
   end
end
var_res = var_total - var_res;
%t = t_best;
%q1 = sum(P(1:(t-1)));
%q2 = sum(P(t:K));
%mu1 = sum(P(1:(t-1)).*(1:(t-1)))/q1;
%var1 = sum(P(1:(t-1)).*(((1:(t-1))-mu1).^2))/q1;
%mu2 = sum(P(t:K).*(t:K))/q2;
%var2 = sum(P(t:K).*(((t:K)-mu2).^2))/q2;
%var_res = q1*var1+q2*var2;
%keyboard
return;
 
