function t = isodata(image)
bins = 256;     %assume 256 colors

h = hist(image(:), bins);
c = [1:bins] / bins;

[m1, i1] = min(h);
[m2, i2] = max(h);
t = round(0.5 * (i1 + i2));

for a = 1:100
    m1 = sum(h(1:t) .* c(1:t)) ./ sum(h(1:t));
    m2 = sum(h(t:end) .* c(t:end)) ./ sum(h(t:end));
    nt = round(bins * 0.5 * (m1 + m2));
    if t == nt
        t=t/bins;
        return;
    end
    t = nt;
end
t=t/bins;