function [T] = otsu_th(I, N)
   [hist_vect, bin_edges] = histcounts(I(:), N);
   
   hist_vect_norm = hist_vect ./ sum(hist_vect,2);
   
   T = zeros(size(I));
    
   for i = 1:size(hist_vect, 2)
       omega_1_k(i) = sum(sum(hist_vect_norm(1:i)));
       mu_1_k(i) = (1/omega_1_k(i)) * sum(sum(i*hist_vect_norm(1:i)));           
              
       omega_2_k(i) = sum(sum(hist_vect_norm(i+1:size(hist_vect,2))));
       mu_2_k(i) = (1/omega_2_k(i)) * sum(sum((i+1)*hist_vect_norm(i+1:size(hist_vect,2))));
       
       sigma_b_k(i) = omega_1_k(i) * omega_2_k(i) * (mu_1_k(i) - mu_2_k(i))^2;
   end
   
  [~, arg_max_sigma_b] = max(sigma_b_k);
  
  th = bin_edges(arg_max_sigma_b + 1);
  
  T = I(:,:,1) > th;
   
end