close all; clear all; clc;

Im=imread('pout.tif');
figure, imshow(Im)
title('Init image')
[nr, nc]=size(Im);
Imsize=nr*nc;
figure
[P,x]=imhist(Im);
bar(x,P);
title('histogram')
p=P/Imsize;

for k=1:256
       c0(1:k)=p(1:k);
       c1(k+1:256)=p(k+1:256);
       W0=sum(c0);
       if W0>0
           W1=1-W0;
           if W1>0
               m0=0;
               for i=1:k
                   m0=m0+i*p(i);
               end
               m0=m0/W0;
               m1=0;
               for i=k+1:256
                   m1=m1+i*p(i);
               end
               m1=m1/W1;
               sb(k)=W0*W1*(m1-m0)*(m1-m0);
           end
       end
end

[peak,Threshold]=min(sb);
Threshold;
for i=0:255
    x(i+1)=i;
end
x1(1)=0;
x1(2)=Threshold;
y1(2)=peak;
figure
plot(x,sb,x1,y1)
title('sb^2')
text(Threshold-50,450,'Threshold position');
Im=im2double(Im);
bw=im2bw(Im,Threshold/255.0);
figure
imshow(bw)