function main

    img=imread(\'lena.jpg\');
    imshow(img);
    img=double(img);
    [m n]=size(img);
    
    Hist=zeros(1,256);
    for i=1:m
        for j=1:n
            Hist(img(i,j)+1)=Hist(img(i,j)+1)+1;    %求直方图
        end
    end

    p=Hist/(m*n);       %直方图概率分布
    uT=sum((1:256).*p(1:256));  %图像亮度均值，其实比真正的均值要大1，所以后面减了1

    sigma_2=zeros(1,256);
    for k=1:256
        sigma_2(k)=(uT*w(k,p)-u(k,p))^2/(w(k,p)*(1-w(k,p)));    %类间方差 
    end
    [tmp index]=max(sigma_2);       %求最大类间方差的索引
    index=index-1;                  %这里索引是1-256，实际图像灰度是0-255，所以减1
    
    imgn=img>index;
    figure;
    imshow(imgn);
    
    function re=w(k,p)      %直方图前k个亮度级的0阶累积矩
        re=sum(p(1:k));   
    end

    function re=u(k,p)      %直方图前k个亮度级的1阶累积矩
        re=sum((1:k).*p(1:k));
    end
        
end