% function thresh_md(a)
tic;
a=imread('kids.tif');
% a=rgb2gray(A);
%  該函數實現最大方差法計算分割門限。
%  輸入參數為灰度圖像，輸出為灰度分割門限
count1=imhist(a);
%  返回圖像矩陣a各個灰度等級象素個數。
[m,n]=size(a);
gs=find(a==0);
mm=size(gs)
N=m*n-sum(sum(find(a==0),1));

L=256;
%  指定圖像灰度等級為256級。
%  d=std2((double(a)+1)/N)
%  err=0.003;
count=count1/N;
%  計算出個灰度出現的概率。
for i=2:L
    if count(i)~=0
        st=i-1;
        break
    end
end
%   找出概率不為0的最小灰度。
for i=L:-1:1
    if count(i)~=0;
        nd=i-1;
        break;
    end
end
%  找出概率不為0的最大灰度。
f=count(st+1:nd+1);
p=st;q=nd-st;
%   p和q分別為灰度起始和結束值。
u=0;
for i=1:q
    u=u+f(i)*(p+i-1);
    ua(i)=u;
end
%  計算圖像的平均灰度。
for i=1:q
    w(i)=sum(f(1:i));
end
%  計算出選擇不同k值時，A區域的概率。
d=(u*w-ua).^2./(w.*(1-w));
%  求出不同k值時類間的方差
[y,tp]=max(d);
%  求出最大方差對應的灰度級
th=tp+p;
toc