function [ th ] = otsu(im)
% Selects a threshold for image usinq Otsu's method
% Usage: [threshold]= otsu(image)
%
[h, bins] = imhist(im);
h = h/sum(h) ;               % normalized histogram, with numbers between 0 and 1 in each bin
L=length(h);

kv = 1:L;
omega = cumsum(h);          % 0-order cumulative moment = omega(k) for all ks
mu_k = cumsum(h .* kv');  % 1st-order cum moment = u(k) for all ks
                                % position k = sum(norm_h(i)*bin(i)) from 1 to k
                                % i.e. position 3 = sum(norm_h(i)*bin(i)) from 1 to 3
mu_T = mu_k(end);

max_sigma_sq = 0;
th = 0 ; 

% sequential calculation
for k = kv
   new_sigma_sq = 0;    
   if (omega(k) ~= 0) 
       if (omega(k) ~= 1)  
            new_sigma_sq = ((mu_T*omega(k) - mu_k(k))^2)/( omega(k)*(1 - omega(k) ) ); % between-class variance
        end
    end
   if (new_sigma_sq > max_sigma_sq)
		max_sigma_sq = new_sigma_sq;
		idx = k;
   end
end

% vector calc
sigma_b_squared = (mu_T * omega - mu_k).^2 ./ (omega .* (1 - omega));
maxval = max(sigma_b_squared);
% find the mean in case there are several maxima
idx = mean(find(sigma_b_squared == maxval));
th = idx - 1;



