%讀入圖片，并灰度化
I=imread（‘鋼筋根數(shù)統(tǒng)計.jpg‘）;
figureimshow（I）;
I=rgb2gray（I）;
%計算圖片的像素點個數(shù)
[mn]=size（I）;
%存放各灰度級出現(xiàn)的次數(shù)
num=zeros（1256）;?
%存放各灰度級的比率
p=zeros（1256）;?
%雙精度化
I=double（I）;?
for?i=1:m?
????for?j=1:n?
???????num（I（ij）+1）=num（I（ij）+1）+1;%統(tǒng)計各灰度級的像素點個數(shù)
????end
end
for?i=1:256?
????p（i）=num（i）/（m*n）;%計算各灰度級出現(xiàn)的比率?
end?
for?i=2:256?
????if?p（i）~=0?
????????st=i+1;%實現(xiàn)尋找出現(xiàn)比率不為0的最小灰度值?
????????break?
????end?
end?
for?i=256:-1:1
????if?p（i）~=0;
????????nd=i-1;%實現(xiàn)找出出現(xiàn)比率不為0的最大灰度值
????????break?
????end?
end

%進行otsu閾值分割
%以下程序?qū)崿F(xiàn)利用最小方差和法找出門閾值?
w=inf;th=0;?
for?t=st:nd%最小非零比率灰度值到最大非零比率灰度值?
????qt1=0;qt2=0;%前景后景像素點比率?
????u1=0;u2=0;%前景后景均值?
????v1=0;v2=0;
????for?i=1:t?
????????qt1=qt1+p（i）;?
????end?
????for?i=1:t?
????????u1=u1+i*p（i）/qt1;?
????end?
????for?i=1:t?
????????v1=v1+（（i-u1）^2）*p（i）/qt1;?
????end?
????for?i=t+1:256?
????????qt2=qt2+p（i）;?
????end?
????for?i=t+1:256?
????????u2=u2+i*p（i）/qt2;?
????end?
????for?i=t+1:256?
????????v2=v2+（（i-u2）^2）*p（i）/qt2;?