大津の二値化法

大津の二値化でしきい値処理した画像の一例

元の画像

大津の二値化法（おおつのにちかほう、Otsu's method、大津の方法などとも呼ばれる）とは、コンピュータビジョンや画像処理において、自動画像しきい値処理を行うために使用される手法。その名はテンプレート:仮リンクにちなむ。最も単純な形式では、アルゴリズムはピクセルをforegroundとbackgroundの2つのクラスに分ける1つの強度しきい値を返す。このしきい値はクラス内の強度分散を最小化することにより、または同等にクラス間の分散を最大化することにより決定される。

大津の二値化はフィッシャーの判別分析の1次元の離散に相当するものであり、テンプレート:仮リンクに関連しており、強度ヒストグラムで行われる大域最適なK平均法と同等である。マルチレベルのしきい値処理へ拡大することは最初の論文で説明されており、以降計算効率の高い実装が提案されている。

大津の二値化

大津の二値化を視覚化したもの

アルゴリズムは、2つのクラスの分散の加重和として定義されるクラス内分散を最小化するしきい値を探す。

${\sigma }_w^2(t)={\omega }_0(t){\sigma }_0^2(t)+{\omega }_1(t){\sigma }_1^2(t)$

加重${\omega }_0$と${\omega }_1$はしきい値$t$により分けられた2つのクラスの確率であり、${\sigma }_0^2$と${\sigma }_1^2$はこれら2つのクラスの分散である。

クラスの確率${\omega }_{0,1}(t)$はヒストグラムの$L$個のビンから計算される。

$\begin{array}{cc}{\omega }_0(t)& =\sum _{i=0}^{t-1}p(i)\\ {\omega }_1(t)& =\sum _{i=t}^{L-1}p(i)\end{array}$

2つのクラスでは、クラス内分散を最小化することはクラス間分散を最大化することと同じである。

$\begin{array}{cc}{\sigma }_b^2(t)& ={\sigma }^2-{\sigma }_w^2(t)={\omega }_0(t)({\mu }_0-{\mu }_T{)}^2+{\omega }_1(t)({\mu }_1-{\mu }_T{)}^2\\ & ={\omega }_0(t){\omega }_1(t){[{\mu }_0(t)-{\mu }_1(t)]}^2\end{array}$

これはクラス確率$\omega$とクラス平均$\mu$で表され、クラス平均${\mu }_0(t)$、${\mu }_1(t)$および${\mu }_T$は、

$\begin{array}{cc}{\mu }_0(t)& =\frac{\sum _{i=0}^{t-1}ip(i)}{{\omega }_0(t)}\\ {\mu }_1(t)& =\frac{\sum _{i=t}^{L-1}ip(i)}{{\omega }_1(t)}\\ {\mu }_T& =\sum _{i=0}^{L-1}ip(i)\end{array}$

である。以上の関係は以下のように簡単に確かめることができる。

$\begin{array}{cc}{\omega }_0{\mu }_0+{\omega }_1{\mu }_1& ={\mu }_T\\ {\omega }_0+{\omega }_1& =1\end{array}$

クラス確率とクラス平均は、繰り返し計算することができる。このアイデアは効果的なアルゴリズムを生む。

アルゴリズム

1. 各強度レベルのヒストグラムと確率を計算する。
2. ${\omega }_i(0)$と${\mu }_i(0)$の初期値を設定する。
3. $t=1,\dots$ 最大強度までの全ての考えうるしきい値で次のステップを行う。
   1. ${\omega }_i$と${\mu }_i$を更新する。
   2. ${\sigma }_b^2(t)$を計算する。
4. 望ましいしきい値は${\sigma }_b^2(t)$が最大となる値である。

MATLABまたはOctaveでの実装

histogramCountsはさまざまなグレーレベルのグレースケール画像（通常は8ビット画像）の256画素のヒストグラムである。levelは画像のしきい値 (double) である。

function level = otsu(histogramCounts)
total = sum(histogramCounts); % total number of pixels in the image 
%% OTSU automatic thresholding
top = 256;
sumB = 0;
wB = 0;
maximum = 0.0;
sum1 = dot(0:top-1, histogramCounts);
for ii = 1:top
    wF = total - wB;
    if wB > 0 && wF > 0
        mF = (sum1 - sumB) / wF;
        val = wB * wF * ((sumB / wB) - mF) * ((sumB / wB) - mF);
        if ( val >= maximum )
            level = ii;
            maximum = val;
        end
    end
    wB = wB + histogramCounts(ii);
    sumB = sumB + (ii-1) * histogramCounts(ii);
end
end

MatlabにはImage Processing Toolboxに組み込み関数graythresh()とmultithresh()があり、それぞれ大津の手法とマルチ大津の手法(Multi Otsu's method)で実装されている。

Pythonでの実装

この実装にはNumPyライブラリが必要である

def compute_otsu_criteria(im, th):
    # create the thresholded image
    thresholded_im = np.zeros(im.shape)
    thresholded_im[im >= th] = 1

    # compute weights
    nb_pixels = im.size
    nb_pixels1 = np.count_nonzero(th)
    weight1 = nb_pixels1 / nb_pixels
    weight0 = 1 - weight1

    # if one the classes is empty, eg all pixels are below or above the threshold, that threshold will not be considered
    # in the search for the best threshold
    if weight1 == 0 or weight0 == 0:
        return np.nan

    # find all pixels belonging to each class
    val_pixels1 = im[thresholded_im == 1]
    val_pixels0 = im[thresholded_im == 0]

    # compute variance of these classes
    var0 = np.var(val_pixels0) if len(val_pixels0) > 0 else 0
    var1 = np.var(val_pixels1) if len(val_pixels1) > 0 else 0

    return weight0 * var0 + weight1 * var1

im = # load your image as a numpy array.
# For testing purposes, one can use for example im = np.random.randint(0,255, size = (50,50))

# testing all thresholds from 0 to the maximum of the image
threshold_range = range(np.max(im)+1)
criterias = [compute_otsu_criteria(im, th) for th in threshold_range]

# best threshold is the one minimizing the Otsu criteria
best_threshold = threshold_range[np.argmin(criterias)]

OpenCVやScikit-imageなどの画像処理専用のPythonライブラリは、アルゴリズムの組み込み実装を提案する。

限界とバリエーション

大津の方法は、ヒストグラムが2つのピークの間に深く鋭い谷を持つ二峰性分布を有するときにうまく機能する。

他の全ての大域的なしきい値処理と同様に、ノイズが大きく、対象のサイズが小さく、明暗が不均一で、クラス内の分散がクラス間の分散よりも大きい場合にパフォーマンスが低下する。このような場合、大津の方法を局所的に適用することが開発されている。

さらに、大津の方法の数学的根拠は、画像のヒストグラムを分散が等しくサイズが等しい2つの正規分布の混合としてモデル化する。しかし、大津のしきい値処理はこれらの仮定を満たさない場合でも満足のいく結果をもたらす可能性がある。同様に統計検定（大津の方法が密接に関連する）は、動作の仮定が完全に満たされていない場合でも正しく実行される。しかし、これらの仮定から生じる深刻な逸脱を無くすために、Kittler-Illingworthの方法などの大津の方法のいくつかのバリエーションが提案されている。

ノイズの多い画像に対するバリエーション

大津の方法の限界に対処するためにさまざまな拡張が開発されてきた。人気のある拡張の1つは、ノイズの多い画像のオブジェクトセグメンテーションのタスクに適した2次元大津の方法である。ここでは、特定のピクセルの強度値をそのすぐ近傍の平均強度と比較することでセグメンテーション結果を改良する。

各ピクセルで、近傍の平均グレーレベル値が計算される。与えられたピクセルのグレーレベルを$L$個の離散値に分割し、平均グレーレベルも同じ$L$個の値に分割する。その後、ピクセルのグレーレベルと近傍の平均のペア$(i,j)$が形成される。各ペアは$L\times L$個の考えうる2次元ビンの1つに属する。ペア$(i,j)$の出現の総数（頻度）$f_{ij}$を画像のピクセルの総数$N$で割ったものは、2次元ヒストグラムで同時確率密度関数を定義する。

$P_{ij}=\frac{f_{ij}}{N},\sum _{i=0}^{L-1}\sum _{j=0}^{L-1}{P}_{ij}=1$

そして、2次元大津の方法は、2次元ヒストグラムに基づいて次のように開発されている。

2つのクラスの確率は、次のように表せる。

$\begin{array}{cc}{\omega }_0& =\sum _{i=0}^{s-1}\sum _{j=0}^{t-1}{P}_{ij}\\ {\omega }_1& =\sum _{i=s}^{L-1}\sum _{j=t}^{L-1}{P}_{ij}\end{array}$

2つのクラスの強度平均ベクトルと合計平均ベクトルは次のように表せる。

$\begin{array}{cc}{\mu }_0& =[{\mu }_{0i},{\mu }_{0j}{]}^T={[\sum _{i=0}^{s-1}\sum _{j=0}^{t-1}i\frac{P_{ij}}{{\omega }_0},\sum _{i=0}^{s-1}\sum _{j=0}^{t-1}j\frac{P_{ij}}{{\omega }_0}]}^T\\ {\mu }_1& =[{\mu }_{1i},{\mu }_{1j}{]}^T={[\sum _{i=s}^{L-1}\sum _{j=t}^{L-1}i\frac{P_{ij}}{{\omega }_1},\sum _{i=s}^{L-1}\sum _{j=t}^{L-1}j\frac{P_{ij}}{{\omega }_1}]}^T\\ {\mu }_T& =[{\mu }_{Ti},{\mu }_{Tj}{]}^T={[\sum _{i=0}^{L-1}\sum _{j=0}^{L-1}i{P}_{ij},\sum _{i=0}^{L-1}\sum _{j=0}^{L-1}j{P}_{ij}]}^T\end{array}$

ほとんどの場合、非対角の確率はわずかであるため、次のことを簡単に確認することができる。

${\omega }_0+{\omega }_1\cong 1$

${\omega }_0{\mu }_0+{\omega }_1{\mu }_1\cong {\mu }_T$

クラス間離散行列は次のように定義される。

$S_b=\sum _{k=0}^1{\omega }_k[({\mu }_k-{\mu }_T)({\mu }_k-{\mu }_T{)}^T]$

離散行列のトレースは次のように表される。

$\begin{array}{cc}& \mathrm{tr}(S_b)\\ =& {\omega }_0[({\mu }_{0i}-{\mu }_{Ti}{)}^2+({\mu }_{0j}-{\mu }_{Tj}{)}^2]+{\omega }_1[({\mu }_{1i}-{\mu }_{Ti}{)}^2+({\mu }_{1j}-{\mu }_{Tj}{)}^2]\\ =& \frac{({\mu }_{Ti}{\omega }_0-{\mu }_i{)}^2+({\mu }_{Tj}{\omega }_0-{\mu }_j{)}^2}{{\omega }_0(1-{\omega }_0)}\end{array}$

ここで

${\mu }_i=\sum _{i=0}^{s-1}\sum _{j=0}^{t-1}i{P}_{ij}$

${\mu }_j=\sum _{i=0}^{s-1}\sum _{j=0}^{t-1}j{P}_{ij}$

1次元の大津の方法と同様に、最適なしきい値$(s,t)$は、$\mathrm{tr}(S_b)$を最大化することで取得される。

アルゴリズム

$s$と$t$は、1次元の大津の方法と同様に繰り返し取得される。$s$と$t$の値は、$\mathrm{tr}(S_b)$の最大値を取得するまで変更される。つまり、

max,s,t = 0;
for ss: 0 to L-1 do
    for tt: 0 to L-1 do
        evaluate tr(S_b);
        if tr(S_b) > max
            max = tr(S,b);
            s = ss;
            t = tt;
        end if
    end for
end for
return s,t;

$\mathrm{tr}(S_b)$を評価するために、高速再起動的プログラミングアルゴリズムを使用して時間パフォーマンスを向上させることができることに留意してください。しかし、動的プログラミングのアプローチを使用しても、2次元大津の方法は依然として時間計算量が非常に複雑である。したがって、計算コストを削減するために多くの研究が行われてきた。

合計領域テーブルを使用して3つのテーブルを作製する場合は、$P_{ij}$を合計し、$i*P_{ij}$を合計し、$j*P_{ij}$を合計し、その時のランタイムの複雑性は(O(N_pixels), O(N_bins*N_bins))の最大値である。しきい値に関して粗い解像度のみが必要な場合はN_binsを減らすことができることに注意してください。

en:Summed-area table参照

Matlabでの実装

関数の入力と出力

histsは、グレースケール値と近傍平均グレースケール値のペアの$256\times 256$の2次元ヒストグラムである。

totalは、所与の画像のペアの数である。これは各方向の2次元ヒストグラムのビンの数により決まる。

thresholdは取得されたしきい値である。

function threshold = otsu_2D(hists, total)
maximum = 0.0;
threshold = 0;
helperVec = 0:255;
mu_t0 = sum(sum(repmat(helperVec',1,256).*hists));
mu_t1 = sum(sum(repmat(helperVec,256,1).*hists));
p_0 = zeros(256);
mu_i = p_0;
mu_j = p_0;
for ii = 1:256
    for jj = 1:256
        if jj == 1
            if ii == 1
                p_0(1,1) = hists(1,1);
            else
                p_0(ii,1) = p_0(ii-1,1) + hists(ii,1);
                mu_i(ii,1) = mu_i(ii-1,1)+(ii-1)*hists(ii,1);
                mu_j(ii,1) = mu_j(ii-1,1);
            end
        else
            p_0(ii,jj) = p_0(ii,jj-1)+p_0(ii-1,jj)-p_0(ii-1,jj-1)+hists(ii,jj); % THERE IS A BUG HERE. INDICES IN MATLAB MUST BE HIGHER THAN 0. ii-1 is not valid
            mu_i(ii,jj) = mu_i(ii,jj-1)+mu_i(ii-1,jj)-mu_i(ii-1,jj-1)+(ii-1)*hists(ii,jj);
            mu_j(ii,jj) = mu_j(ii,jj-1)+mu_j(ii-1,jj)-mu_j(ii-1,jj-1)+(jj-1)*hists(ii,jj);
        end

        if (p_0(ii,jj) == 0)
            continue;
        end
        if (p_0(ii,jj) == total)
            break;
        end
        tr = ((mu_i(ii,jj)-p_0(ii,jj)*mu_t0)^2 + (mu_j(ii,jj)-p_0(ii,jj)*mu_t1)^2)/(p_0(ii,jj)*(1-p_0(ii,jj)));

        if ( tr >= maximum )
            threshold = ii;
            maximum = tr;
        end
    end
end
end

不平衡な画像に対するバリエーション

画像のクラスのグレーレベルは正規分布とみなすことができるが、サイズや分散が等しくない場合、大津のアルゴリズムの仮定は満たされない。Kittler-Illingworthのアルゴリズム（最小エラーしきい値処理とも呼ばれる）はこのような場合を処理するための大津の方法のバリエーションである。このアルゴリズムを数学的に説明する方法はいくつかある。それらの1つは、試験される各しきい値について結果の2値画像の正規分布のパラメータがデータが与えられた最尤推定により推定されたことを考慮することである。

このアルゴリズムは大津の方法よりも優れているように見えることがあるが、推定される新しいパラメータが導入されるためアルゴリズムが過剰パラメータ化されて不安定になる可能性がある。大津の方法からの仮定が少なくとも部分的に有効であると考えられる多くの場合において、オッカムの剃刀に従ってKittler-Illingworthのアルゴリズムよりも大津の方法を優先する方が好ましい場合がある。

外部リンク

• Implementation of Otsu's thresholding method as GIMP-plugin using Script-Fu (a Scheme-based language)
• 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
• Otsu Thresholding with scikit-image in Python