大津算法
![本页使用了标题或全文手工转换](http://images.weserv.nl/?url=//upload.wikimedia.org/wikipedia/commons/thumb/c/cd/Zh_conversion_icon_m.svg/35px-Zh_conversion_icon_m.svg.png)
![](http://images.weserv.nl/?url=//upload.wikimedia.org/wikipedia/commons/thumb/4/46/Image_processing_post_otsus_algorithm.jpg/220px-Image_processing_post_otsus_algorithm.jpg)
![](http://images.weserv.nl/?url=//upload.wikimedia.org/wikipedia/commons/thumb/4/4b/Image_processing_pre_otsus_algorithm.jpg/220px-Image_processing_pre_otsus_algorithm.jpg)
在计算机视觉和图像处理领域,大津二值化法(Otsu's method)是由大津展之命名的一种用于自动图像阈值分割的方法。[1] 该方法用于对基于聚类的图像进行二值化,即将一个灰度图像转化为二值图像。算法假定图像根据双模直方图(前景像素和背景像素)包含两类像素,并计算出能将这两类像素分开的最佳阈值,使得类内方差最小,或等效地,类间方差最大。[2] 因此,大津二值化法是一种一维离散的費舍爾判别分析的类似方法,与Jenks优化方法相关,并且等同于在强度直方图上执行的全局最优k-means算法。 原始论文中还描述了多级阈值扩展方法,[3] 并且后来提出了计算效率高的实现方法。[4][5]
方法
在大津算法中,我们穷举搜索能使类内方差最小的阈值,定义为两个类的方差的加权和:
权重 是被阈值 分开的两个类的概率,而 是这两个类的方差。
大津证明了最小化类内方差和最大化类间方差是相同的:[2]
用类概率 和类均值 来表示。
类概率 用阈值为 的直方图计算:
而类均值 为:
其中 为第 个直方图面元中心的值。 同样的,你可以对大于 的面元求出右侧直方图的 与 。
类概率和类均值可以迭代计算。这个想法会产生一个有效的算法。
大津算法得出了0:1范围上的一个阈值。这个阈值用于图像中出现的像素强度的动态范围。例如,若图像只包含155到255之间的像素强度,大津阈值0.75会映射到灰度阈值230(而不是192,因为图像包含的像素不是0–255全范围的)。常见的摄影图像往往包含全范围的像素强度,就会让这一点有争论,但其他应用会对这点区别很敏感。[6]
算法
- 计算每个强度级的直方图和概率
- 设置 和 的初始值
- 遍历所有可能的阈值 最大强度
- 更新 和
- 计算
- 所需的阈值对应于最大的
- 你可以计算两个最大值(和两个对应的)。 是最大值而 是更大的或相等的最大值
- 所需的阈值 =
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)可以很容易地转变成多线程的方法。
参考文献
- ^ 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.0 2.1 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.
- ^ 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.
- ^ Liao, Ping-Sung. A fast algorithm for multilevel thresholding (PDF). J. Inf. Sci. Eng. 2001, 17 (5): 713–727. S2CID 9609430. doi:10.6688/JISE.2001.17.5.1. (原始内容 (PDF)存档于2019-06-24).
- ^ Huang, Deng-Yuan. Optimal multi-level thresholding using a two-stage Otsu optimization approach. Pattern Recognition Letters. 2009, 30 (3): 275–284. Bibcode:2009PaReL..30..275H. doi:10.1016/j.patrec.2008.10.003.
- ^ 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(页面存档备份,存于互联网档案馆)