大津演算法

在電腦視覺和圖像處理領域，大津二值化法（Otsu's method）是由大津展之（英語：Nobuyuki Otsu）命名的一種用於自動圖像閾值分割的方法。^[1] 該方法用於對基於聚類的圖像進行二值化，即將一個灰度圖像轉化為二值圖像。演算法假定圖像根據雙模直方圖（前景像素和背景像素）包含兩類像素，並計算出能將這兩類像素分開的最佳閾值，使得類內方差最小，或等效地，類間方差最大。^[2] 因此，大津二值化法是一種一維離散的費舍爾判別分析的類似方法，與Jenks最佳化方法（英語：Jenks natural breaks optimization）相關，並且等同於在強度直方圖上執行的全域最佳k-means演算法。原始論文中還描述了多級閾值擴充方法，^[3] 並且後來提出了計算效率高的實現方法。^[4]^[5]

方法

在大津演算法中，我們窮舉搜尋能使類內方差最小的閾值，定義為兩個類的方差的加權和：

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

權重 $\omega _{i}$ 是被閾值 $t$ 分開的兩個類的概率，而 $\sigma _{i}^{2}$ 是這兩個類的方差。

大津證明了最小化類內方差和最大化類間方差是相同的：^[2]

\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}

用類概率 $\omega _{i}$ 和類均值 $\mu _{i}$ 來表示。

類概率 $\omega _{1}(t)$ 用閾值為 $t$ 的直方圖計算：

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

而類均值 $\mu _{1}(t)$ 為：

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

其中 $x(i)$ 為第 $i$ 個直方圖面元中心的值。同樣的，你可以對大於 $t$ 的面元求出右側直方圖的 $\omega _{2}(t)$ 與 $\mu _{2}$ 。

類概率和類均值可以迭代計算。這個想法會產生一個有效的演算法。

大津演算法得出了0:1範圍上的一個閾值。這個閾值用於圖像中出現的像素強度的動態範圍。例如，若圖像只包含155到255之間的像素強度，大津閾值0.75會對映到灰度閾值230（而不是192，因為圖像包含的像素不是0–255全範圍的）。常見的攝影圖像往往包含全範圍的像素強度，就會讓這一點有爭論，但其他應用會對這點區別很敏感。^[6]

演算法

計算每個強度級的直方圖和概率
設置 $\omega _{i}(0)$ 和 $\mu _{i}(0)$ 的初始值
遍歷所有可能的閾值 $t=1\ldots$ $t=1\ldots$ 最大強度
1. 更新 $\omega _{i}$ 和 $\mu _{i}$
2. 計算 $\sigma _{b}^{2}(t)$
所需的閾值對應於最大的 $\sigma _{b}^{2}(t)$
你可以計算兩個最大值（和兩個對應的）。 $\sigma _{b1}^{2}(t)$ 是最大值而 $\sigma _{b2}^{2}(t)$ 是更大的或相等的最大值
所需的閾值 = ${\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）可以很容易地轉變成多線程的方法。

參考文獻

^ 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（頁面存檔備份，存於互聯網檔案館）

[Mehmet-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.

[Otsu-2] 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.

[multi-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] 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）.

[5] 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.

[6] Q&A Forum. Stack Overflow. [2015-10-20]. （原始內容存檔於2015-11-14）.

[1]

[2]

[3]

[4]

[5]

[6]