大津算法

方法

${\displaystyle \sigma _{w}^{2}(t)=\omega _{1}(t)\sigma _{1}^{2}(t)+\omega _{2}(t)\sigma _{2}^{2}(t)}$

${\displaystyle \sigma _{b}^{2}(t)=\sigma ^{2}-\sigma _{w}^{2}(t)=\omega _{1}(t)\omega _{2}(t)\left[\mu _{1}(t)-\mu _{2}(t)\right]^{2}}$

${\displaystyle \omega _{1}(t)=\Sigma _{0}^{t}p(i)}$

${\displaystyle \mu _{1}(t)=\left[\Sigma _{0}^{t}p(i)\,x(i)\right]/\omega _{1}}$

算法

1. 计算每个强度级的直方图和概率
2. 设置 ${\displaystyle \omega _{i}(0)}$ ${\displaystyle \mu _{i}(0)}$  的初始值
3. 遍历所有可能的阈值 ${\displaystyle t=1\ldots }$  最大强度
1. 更新 ${\displaystyle \omega _{i}}$ ${\displaystyle \mu _{i}}$
2. 计算 ${\displaystyle \sigma _{b}^{2}(t)}$
4. 所需的阈值对应于最大的 ${\displaystyle \sigma _{b}^{2}(t)}$
5. 你可以计算两个最大值（和两个对应的）。${\displaystyle \sigma _{b1}^{2}(t)}$  是最大值而 ${\displaystyle \sigma _{b2}^{2}(t)}$  是更大的或相等的最大值
6. 所需的阈值 = ${\displaystyle {\frac {{\text{threshold}}_{1}+{\text{threshold}}_{2}}{2}}}$

JavaScript实现

function otsu(histogram, total) {
var sum = 0;
for (var i = 1; i < 256; ++i)
sum += i * histogram[i];
var sumB = 0;
var wB = 0;
var wF = 0;
var mB;
var mF;
var max = 0.0;
var between = 0.0;
var threshold1 = 0.0;
var threshold2 = 0.0;
for (var i = 0; i < 256; ++i) {
wB += histogram[i];
if (wB == 0)
continue;
wF = total - wB;
if (wF == 0)
break;
sumB += i * histogram[i];
mB = sumB / wB;
mF = (sum - sumB) / wF;
between = wB * wF * (mB - mF) * (mB - mF);
if ( between >= max ) {
threshold1 = i;
if ( between > max ) {
threshold2 = i;
}
max = between;
}
}
return ( threshold1 + threshold2 ) / 2.0;
}


C实现

 1 unsigned char otsu_threshold( int* histogram, int pixel_total )
2 {
3     unsigned int sumB = 0;
4     unsigned int sum1 = 0;
5     float wB = 0.0f;
6     float wF = 0.0f;
7     float mF = 0.0f;
8     float max_var = 0.0f;
9     float inter_var = 0.0f;
10     unsigned char threshold = 0;
11     unsigned short index_histo = 0;
12
13     for ( index_histo = 1; index_histo < 256; ++index_histo )
14     {
15         sum1 += index_histo * histogram[ index_histo ];
16     }
17
18     for (index_histo = 1; index_histo < 256; ++index_histo)
19     {
20         wB = wB + histogram[ index_histo ];
21         wF = pixel_total - wB;
22         if ( wB == 0 || wF == 0 )
23         {
24             continue;
25         }
26         sumB = sumB + index_histo * histogram[ index_histo ];
27         mF = ( sum1 - sumB ) / wF;
28         inter_var = wB * wF * ( ( sumB / wB ) - mF ) * ( ( sumB / wB ) - mF );
29         if ( inter_var >= max_var )
30         {
31             threshold = index_histo;
32             max_var = inter_var;
33         }
34     }
35
36     return threshold;
37 }


MATLAB实现

total为给定图像中的像素数。 histogramCounts为灰度图像不同灰度级（典型的8位图像）的256元直方图。 level为图像的阈值（double）。

function level = otsu(histogramCounts, total)
%% OTSU automatic thresholding method
sumB = 0;
wB = 0;
maximum = 0.0;
threshold1 = 0.0;
threshold2 = 0.0;
sum1 = sum((1:256).*histogramCounts.'); % the above code is replace with this single line
for ii=1:256
wB = wB + histogramCounts(ii);
if (wB == 0)
continue;
end
wF = total - wB;
if (wF == 0)
break;
end
sumB = sumB +  ii * histogramCounts(ii);
mB = sumB / wB;
mF = (sum1 - sumB) / wF;
between = wB * wF * (mB - mF) * (mB - mF);
if ( between >= maximum )
threshold1 = ii;
if ( between > maximum )
threshold2 = ii;
end
maximum = between;
end
end
level = (threshold1 + threshold2 )/(2);
end


function  [threshold_otsu] = Thredsholding_Otsu( Image)
%Intuition:
%(1)pixels are divided into two groups
%(2)pixels within each group are very similar to each other
%   Parameters:
%   t : threshold
%   r : pixel value ranging from 1 to 255
%   q_L, q_H : the number of lower and higher group respectively
%   sigma : group variance
%   miu : group mean
%   Author: Lei Wang
%   Date  : 22/09/2013
%   References : Wikepedia,
%   This is my original work

nbins = 256;
counts = imhist(Image,nbins);
p = counts / sum(counts);

for t = 1 : nbins
q_L = sum(p(1 : t));
q_H = sum(p(t + 1 : end));
miu_L = sum(p(1 : t) .* (1 : t)') / q_L;
miu_H = sum(p(t + 1 : end) .* (t + 1 : nbins)') / q_H;
sigma_b(t) = q_L * q_H * (miu_L - miu_H)^2;
end

[~,threshold_otsu] = max(sigma_b(:));
end


参考文献

1. ^ M. Sezgin and B. Sankur. Survey over image thresholding techniques and quantitative performance evaluation. Journal of Electronic Imaging. 2004, 13 (1): 146–165. doi:10.1117/1.1631315.
2. Nobuyuki Otsu. A threshold selection method from gray-level histograms. IEEE Trans. Sys., Man., Cyber. 1979, 9 (1): 62–66. doi:10.1109/TSMC.1979.4310076.
3. ^ Ping-Sung Liao and Tse-Sheng Chen and Pau-Choo Chung. A Fast Algorithm for Multilevel Thresholding. J. Inf. Sci. Eng. 2001, 17 (5): 713–727.
4. ^ Q&A Forum. Stack Overflow. [2015-10-20]. （原始内容存档于2015-11-14）.