大津算法

在计算机视觉和图像处理中，大津二值化法用来自动对基于聚类的图像进行二值化，^[1] 或者说，将一个灰度图像退化为二值图像。该算法以大津展之（英语：Nobuyuki Otsu）命名。算法假定该图像根据双模直方图（前景像素和背景像素）把包含两类像素，于是它要计算能将两类分开的最佳阈值，使得它们的类内方差最小；由于两两平方距离恒定，所以即它们的类间方差最大。^[2] 因此，大津二值化法粗略的来说就是一维費舍爾判别分析的离散化模拟。

使用大津算法二值化的图像

原图

原始方法的多级阈值扩展称为多大津算法。^[3]

方法

在大津算法中，我们穷举搜索能使类内方差最小的阈值，定义为两个类的方差的加权和：

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

权重 ${\displaystyle \omega _{i}}$  是被阈值 ${\displaystyle t}$  分开的两个类的概率，而 ${\displaystyle \sigma _{i}^{2}}$  是这两个类的方差。

大津证明了最小化类内方差和最大化类间方差是相同的：^[2]

${\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 _{i}}$  和类均值 ${\displaystyle \mu _{i}}$  来表示。

类概率 ${\displaystyle \omega _{1}(t)}$  用阈值为 ${\displaystyle t}$  的直方图计算：

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

而类均值 ${\displaystyle \mu _{1}(t)}$  为：

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

其中 ${\displaystyle x(i)}$  为第 ${\displaystyle i}$  个直方图面元中心的值。 同样的，你可以对大于 ${\displaystyle t}$  的面元求出右侧直方图的 ${\displaystyle \omega _{2}(t)}$  与 ${\displaystyle \mu _{2}}$ 。

类概率和类均值可以迭代计算。这个想法会产生一个有效的算法。

大津算法得出了0:1范围上的一个阈值。这个阈值用于图像中出现的像素强度的动态范围。例如，若图像只包含155到255之间的像素强度，大津阈值0.75会映射到灰度阈值230（而不是192，因为图像包含的像素不是0–255全范围的）。常见的摄影图像往往包含全范围的像素强度，就会让这一点有争论，但其他应用会对这点区别很敏感。^[4]

算法

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实现

注：输入参数total是给定图像中的像素数。输入参数histogram是灰度图像不同灰度级（典型的8位图像）的256元直方图。此函数输出图像的阈值。

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实现

unsigned char otsu_threshold( int* histogram, int pixel_total ) {  unsigned int sumB = 0;  unsigned int sum1 = 0;  float wB = 0.0f;  float wF = 0.0f;  float mF = 0.0f;  float max_var = 0.0f;  float inter_var = 0.0f;  unsigned char threshold = 0;  unsigned short index_histo = 0;   for ( index_histo = 1; index_histo < 256; ++index_histo )  {  sum1 += index_histo * histogram[ index_histo ];  }   for (index_histo = 1; index_histo < 256; ++index_histo)  {  wB = wB + histogram[ index_histo ];  wF = pixel_total - wB;  if ( wB == 0 || wF == 0 )  {  continue;  }  sumB = sumB + index_histo * histogram[ index_histo ];  mF = ( sum1 - sumB ) / wF;  inter_var = wB * wF * ( ( sumB / wB ) - mF ) * ( ( sumB / wB ) - mF );  if ( inter_var >= max_var )  {  threshold = index_histo;  max_var = inter_var;  }  }   return threshold; } 

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 

另一种方法是用向量化方法（可以很容易地转换为便于GPU处理Python矩阵数组版本）

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,  % for multi children Otsu method, please visit : https://drive.google.com/file/d/0BxbR2jt9XyxteF9fZ0NDQ0dKQkU/view?usp=sharing % 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 

该实现有一点点冗余的计算。但由于大津算法快速，这个实现版本是可以接受，并且易于理解的。因为在一些环境中，如果使用向量化的形式，可以更快地运算循环。大津二值化法使用架构最小的堆－孩子标签（heap-children label）可以很容易地转变成多线程的方法。

参考文献

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）. 

外部链接

• Lecture notes on thresholding（页面存档备份，存于互联网档案馆） – covers the Otsu method.
• A plugin for ImageJ（页面存档备份，存于互联网档案馆） using Otsu's method to do the threshold.
• A full explanation of Otsu's method（页面存档备份，存于互联网档案馆） with a working example and Java implementation.
• Implementation of Otsu's method（页面存档备份，存于互联网档案馆） in ITK
• Otsu Thresholding in C#（页面存档备份，存于互联网档案馆） A straightforward C# implementation with explanation.
• Otsu's method using MATLAB（页面存档备份，存于互联网档案馆）