简介

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

方法

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

σ w 2 ( t ) = ω 1 ( t ) σ 1 2 ( t ) + ω 2 ( t ) σ 2 2 ( t ) {displaystyle sigma _{w}^{2}(t)=omega _{1}(t)sigma _{1}^{2}(t)+omega _{2}(t)sigma _{2}^{2}(t)}

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

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

σ b 2 ( t ) = σ 2 − σ w 2 ( t ) = ω 1 ( t ) ω 2 ( t ) 2 {displaystyle sigma _{b}^{2}(t)=sigma ^{2}-sigma _{w}^{2}(t)=omega _{1}(t)omega _{2}(t)left^{2}}

用类概率 ω i {displaystyle omega _{i}} 和类均值 μ i {displaystyle mu _{i}} 来表示。

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

ω 1 ( t ) = Σ 0 t p ( i ) {displaystyle omega _{1}(t)=Sigma _{0}^{t}p(i)}

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

μ 1 ( t ) = / ω 1 {displaystyle mu _{1}(t)=left/omega _{1}}

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

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

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

算法

计算每个强度级的直方图和概率

设置 ω i ( 0 ) {displaystyle omega _{i}(0)} 和 μ i ( 0 ) {displaystyle mu _{i}(0)} 的初始值

遍历所有可能的阈值 t = 1 … {displaystyle t=1ldots } 最大强度

更新 ω i {displaystyle omega _{i}} 和 μ i {displaystyle mu _{i}}

计算 σ b 2 ( t ) {displaystyle sigma _{b}^{2}(t)}

所需的阈值对应于最大的 σ b 2 ( t ) {displaystyle sigma _{b}^{2}(t)}

你可以计算两个最大值（和两个对应的）。 σ b 1 2 ( t ) {displaystyle sigma _{b1}^{2}(t)} 是最大值而 σ b 2 2 ( t ) {displaystyle sigma _{b2}^{2}(t)} 是更大的或相等的最大值

所需的阈值 = threshold 1 + threshold 2 2 {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;    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 = 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;    }    for (index_histo = 1; index_histo = 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 methodsumB = 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 linefor 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;    endendlevel = (threshold1 + threshold2 )/(2);end

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

function   = 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 worknbins = 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 = max(sigma_b(:));end

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