大津の二値化法

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

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

## 大津の二値化

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

:\sigma^2_w(t)=\omega_0(t)\sigma^2_0(t)+\omega_1(t)\sigma^2_1(t)
加重\omega_{0}と\omega_{1}はしきい値tにより分けられた2つのクラスの確率であり、\sigma^2_{0}と\sigma^2_{1}はこれら2つのクラスの分散である。

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

: 
\begin{align}
\omega_0(t) & =\sum_{i=0}^{t-1} p(i)\\[4pt]
\omega_1(t) & =\sum_{i=t}^{L-1} p(i)
\end{align}


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

:
\begin{align}
\sigma^2_b(t) & =\sigma^2-\sigma^2_w(t)=\omega_0(t)(\mu_0-\mu_T)^2+\omega_1(t)(\mu_1-\mu_T)^2 \\
& =\omega_0(t) \omega_1(t) \left[\mu_0(t)-\mu_1(t)\right]^2
\end{align}


これはクラス確率\omegaとクラス平均\muで表され、クラス平均\mu_0(t)、\mu_1(t)および\mu_Tは、

:
\begin{align}
\mu_0(t) & = \frac{\sum_{i=0}^{t-1} i p(i)}{\omega_0(t)} \\[4pt]
\mu_1(t) & = \frac{\sum_{i=t}^{L-1} i p(i)}{\omega_1(t)} \\
\mu_T & = \sum_{i=0}^{L-1} ip(i)
\end{align}


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

:
\begin{align}
\omega_0\mu_0+\omega_1\mu_1 & = \mu_T \\
\omega_0+\omega_1 & =1
\end{align}


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

### アルゴリズム
# 各強度レベルのヒストグラムと確率を計算する。
# \omega_i(0)と\mu_i(0)の初期値を設定する。
# t = 1, \ldots  最大強度までの全ての考えうるしきい値で次のステップを行う。
## \omega_iと\mu_iを更新する。
## \sigma^2_b(t)を計算する。
# 望ましいしきい値は\sigma^2_b(t)が最大となる値である。

### MATLABまたはOctaveでの実装
**<code>histogramCounts</code>**はさまざまなグレーレベルのグレースケール画像（通常は8ビット画像）の256画素のヒストグラムである。`level`は画像のしきい値 (double) である。


```matlab
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ライブラリが必要である。


```python
import numpy as np

def compute_otsu_criteria(im, th):
    # 閾値処理された画像（二値化画像）を作成する
    thresholded_im = np.zeros(im.shape)
    thresholded_im[im >= th] = 1

    # 重みを計算する
    nb_pixels = im.size
    nb_pixels1 = np.count_nonzero(th)
    weight1 = nb_pixels1 / nb_pixels
    weight0 = 1 - weight1

    # クラスの1つが空の場合（例：すべてのピクセルが閾値以下または以上の場合）、
    # その閾値は最適な閾値の探索において考慮されない
    if weight1 == 0 or weight0 == 0:
        return np.nan

    # 各クラスに属するすべてのピクセルを見つける
    val_pixels1 = im[thresholded_im == 1]
    val_pixels0 = im[thresholded_im == 0]

    # これらのクラスの分散を計算する
    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 =  # 画像をnumpy配列として読み込む。
# テスト目的の場合、例えば im = np.random.randint(0,255, size = (50,50)) などを使用できる

# 0から画像の最大値までのすべての閾値をテストする
threshold_range = range(np.max(im) + 1)
criterias = [compute_otsu_criteria(im, th) for th in threshold_range]

# 最適な閾値は、大津の基準（クラス内分散の重み付き和）を最小化するものである
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, \qquad \sum_{i=0}^{L-1}\sum_{j=0}^{L-1} P_{ij}=1

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

2つのクラスの確率は、次のように表せる。
:
\begin{align}
\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{align}


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

:
\begin{align}
\mu_0 & =[\mu_{0i}, \mu_{0j}]^T = \left[\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} \right]^T \\
\mu_1 & =[\mu_{1i}, \mu_{1j}]^T = \left[\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} \right]^T \\
\mu_T & =[\mu_{Ti}, \mu_{Tj}]^T = \left[\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}\right]^T
\end{align}


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

:\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{align}
& \operatorname{tr}(S_b) \\[4pt]
= {} & \omega_0[(\mu_{0i}-\mu_{Ti})^2 + (\mu_{0j}-\mu_{Tj})^2] + \omega_1[(\mu_{1i}-\mu_{Ti})^2 + (\mu_{1j}-\mu_{Tj})^2] \\[4pt]
= {} & \frac{(\mu_{Ti}\omega_0-\mu_i)^2 + (\mu_{Tj}\omega_0-\mu_j)^2}{\omega_0(1-\omega_0)}
\end{align}


ここで

:\mu_i = \sum_{i=0}^{s-1}\sum_{j=0}^{t-1}iP_{ij}
:\mu_j = \sum_{i=0}^{s-1}\sum_{j=0}^{t-1}jP_{ij}

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

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


```matlab
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;
```


\operatorname{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`は取得されたしきい値である。


```matlab
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