%%%  下面是最大方差法計算灰度分割門限的代碼：(代碼以子函數形式給出)
function th = thresh_md(a);
%  該函數實現最大方差法計算分割門限
%  輸入參數為灰度圖像，輸出為灰度分割門限
count = imhist(a);
%  返回圖像矩陣a各個灰度等級像素個數
[m,n] = size(a);
N = m*n - sum(sum(find(a==0),1));
L = 256;
%  計算出現概率不為0的最小灰度
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
f = count(st+1:nd+1);
%  p和q分別為灰度起始和結束值
p = st;
q = nd-st;
% 計算圖像的平均灰度
u = 0;
for i=1:q
u=u+f(i)*(p+i-1);
ua(i)=u;
end
% 計算出選擇不同k值時，A區域的概率
for i=1:q
w(i)=sum(f(1:i));
end
% 求出不同k值時類間的方差
d = (u*w-ua).^2./(w.*(1-w));
% 求出最大方差對應的灰度級
[y,tp] = max(d);
th = tp+p;