function [im2,T] = dootsuthreshold(difim, nbins, mask, usezeros)

%
% Otsu method
%


if (nargin < 2)
    nbins = 256;
end
if (nargin < 3)
    mask=1;
end
if (nargin < 4)
    usezeros=1;
end

data = removezeros(makelinear((difim+1000).*mask))-1000;

[h,hc] = hist(data, nbins);
if (usezeros==0)
    h(1)=0;
end
T = hc(otsuthreshold(h));

im2 = 0+(difim>=T);