function theta = binarize(image)

  % get size of image
  %%%%%%%%%%%%%%%%%%%

  s=size(image);
  height=s(1);
  width=s(2);

  % calculate probabilities for single pixel values
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  % get histogram with 256 bins
  h= hist(image(:),[0:255]);

  % devide each element by sum of all histogram entries
  p_v= h./sum(h);

  % calculate mean value of histogram
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  mu = [0:255]*p_v';

  % calculate p(omega_1) and p(omega_2) for all l=0:255
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  p_o1=tril(ones(256))*p_v';
  p_o2=ones(256,1)-p_o1;

  % calculate mean for both classes
  % => mu_t is a vector for class t with mu_t[i] = mu_t for threshold i 
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
  v_1 = tril(ones(256,1)*p_v)*[0:255]';
  mu_1 = v_1./p_o1; 
  v_2 = triu(ones(256,1)*p_v,1)*[0:255]';
  mu_2 = v_2./p_o2;

  % calculate J
  % => j is a vector of 256 elements with j[i] = J for threshold i 
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

  j = p_o1 .* ((mu_1-mu).*(mu_1-mu)) + p_o2 .* ((mu_2-mu).*(mu_2-mu));
  [val,id]=max(j);
  
  % return threshold
  theta= id-1;
end
