%本程序是利用最大類間方差算法求解自適應閾值，對圖像進行分割
clear;
close all
warning off;

I=imread('11.jpg');
%rgb轉灰度
if isrgb(I)==1
    I_gray=rgb2gray(I);
else
    I_gray=I;
end
figure,imshow(I_gray);
I_double=double(I_gray);%轉化為雙精度
[wid,len]=size(I_gray);
colorlevel=256;  %灰度級
hist=zeros(colorlevel,1);%直方圖
%threshold=128; %初始閾值

%計算直方圖
for i=1:wid
    for j=1:len
        m=I_gray(i,j)+1;
        hist(m)=hist(m)+1;
    end
end
hist=hist/(wid*len);%直方圖歸一化
miuT=0;
for m=1:colorlevel
    miuT=miuT+(m-1)*hist(m);
end
xigmaB2=0;

for mindex=1:colorlevel
    threshold=mindex-1;
    omega1=0;
    omega2=0;
    for m=1:threshold-1
          omega1=omega1+hist(m);
    end
    omega2=1-omega1;
    miu1=0;
    miu2=0;

    for m=1:colorlevel
        if m<threshold
           miu1=miu1+(m-1)*hist(m);
        else
           miu2=miu2+(m-1)*hist(m);
        end
    end
    miu1=miu1/omega1;
    miu2=miu2/omega2;
    xigmaB21=omega1*(miu1-miuT)^2+omega2*(miu2-miuT)^2;
    xigma(mindex)=xigmaB21;
    if xigmaB21>xigmaB2
        finalT=threshold;
        xigmaB2=xigmaB21;
    end
end
%  fT=finalT/255 %閾值歸一化
%  T=graythresh(I_gray)%matlab函數求閾值
for i=1:wid
    for j=1:len
        if I_double(i,j)>finalT
            bin(i,j)=1;
        else
            bin(i,j)=0;
        end
    end
end
figure,imshow(bin);
figure,plot(1:colorlevel,xigma);