大津の二値化法のソースを表示

[[Image:Image processing post otsus algorithm.jpg|thumb|right|大津の二値化でしきい値処理した画像の一例]]
[[Image:Image processing pre otsus algorithm.jpg|thumb|right|元の画像]]
'''大津の二値化法'''（おおつのにちかほう、Otsu's method、'''大津の方法'''などとも呼ばれる）とは、[[コンピュータビジョン]]や[[画像処理]]において、自動画像[[しきい値処理 (画像処理)|しきい値処理]]を行うために使用される手法<ref name=Mehmet>{{cite journal|author1=M. Sezgin |author2=B. Sankur |name-list-style=amp | year = 2004| title = Survey over image thresholding techniques and quantitative performance evaluation| volume = 13| issue = 1| pages = 146-165| journal = Journal of Electronic Imaging| doi = 10.1117/1.1631315}}</ref>。その名は{{仮リンク|大津展之|en|Nobuyuki Otsu}}にちなむ。最も単純な形式では、アルゴリズムはピクセルをforegroundとbackgroundの2つのクラスに分ける1つの強度しきい値を返す。このしきい値はクラス内の強度分散を最小化することにより、または同等にクラス間の分散を最大化することにより決定される<ref name=Otsu>{{cite journal| author = Nobuyuki Otsu| year = 1979| title = A threshold selection method from gray-level histograms| volume = 9| issue = 1| pages = 62-66| journal = IEEE Trans. Sys. Man. Cyber.| doi = 10.1109/TSMC.1979.4310076}}</ref>。

大津の二値化はフィッシャーの判別分析の1次元の離散に相当するものであり、{{仮リンク|ジェンクス最適化法|en|Jenks optimization method}}に関連しており、強度ヒストグラムで行われる大域最適な[[K平均法]]と同等である<ref>{{Cite journal|last=Liu|first=Dongju|date=2009|title=Otsu method and K-means|journal=Ninth International Conference on Hybrid Intelligent Systems IEEE|volume=1|pages=344-349}}</ref>。マルチレベルのしきい値処理へ拡大することは最初の論文で説明されており<ref name="Otsu" />、以降計算効率の高い実装が提案されている<ref>{{Cite journal|last=Liao|first=Ping-Sung|date=2001|title=A fast algorithm for multilevel thresholding|url=https://pdfs.semanticscholar.org/b809/14dbbee9f6b2455742d8117417731e6ecf12.pdf|archiveurl=https://web.archive.org/web/20190624180800/https://pdfs.semanticscholar.org/b809/14dbbee9f6b2455742d8117417731e6ecf12.pdf|url-status=dead|archivedate=2019-06-24|journal=J. Inf. Sci. Eng.|volume=17|issue=5|pages=713-727}}</ref><ref>{{Cite journal|last=Huang|first=Deng-Yuan|date=2009|title=Optimal multi-level thresholding using a two-stage Otsu optimization approach|journal=Pattern Recognition Letters|volume=30|issue=3|pages=275-284|doi=10.1016/j.patrec.2008.10.003}}</ref>。

==大津の二値化==
[[File:Otsu's Method Visualization.gif|thumb|大津の二値化を視覚化したもの]]
アルゴリズムは、2つのクラスの分散の加重和として定義されるクラス内分散を最小化するしきい値を探す。

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

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

: <math>
\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}
</math>

2つのクラスでは、クラス内分散を最小化することはクラス間分散を最大化することと同じである<ref name=Otsu/>。

:<math>
\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}
</math>

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

:<math>
\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}
</math>

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

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

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

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

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


<syntaxhighlight lang="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
</syntaxhighlight>

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

=== Pythonでの実装 ===
この実装には[[NumPy]]ライブラリが必要である<syntaxhighlight lang="python">
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)]
    
</syntaxhighlight>[[OpenCV]]や[[Scikit-image]]などの画像処理専用のPythonライブラリは、アルゴリズムの組み込み実装を提案する。

==限界とバリエーション==
大津の方法は、ヒストグラムが2つのピークの間に深く鋭い谷を持つ二峰性分布を有するときにうまく機能する<ref>{{Cite journal |last=Kittler |first=J. |last2=Illingworth |first2=J. |date=September 1985 |title=On threshold selection using clustering criteria |url=https://doi.org/10.1109/tsmc.1985.6313443 |journal=IEEE Transactions on Systems, Man, and Cybernetics |volume=SMC-15 |issue=5 |pages=652-655 |doi=10.1109/tsmc.1985.6313443 |issn=0018-9472}}</ref>。

他の全ての大域的なしきい値処理と同様に、ノイズが大きく、対象のサイズが小さく、明暗が不均一で、クラス内の分散がクラス間の分散よりも大きい場合にパフォーマンスが低下する<ref name="lee1990comparative">{{cite journal |author=Lee, Sang Uk and Chung, Seok Yoon and Park, Rae Hong |year=1990 |title=A comparative performance study of several global thresholding techniques for segmentation |journal=Computer Vision, Graphics, and Image Processing |volume=52 |issue=2 |pages=171-190 |doi=10.1016/0734-189x(90)90053-x}}</ref>。このような場合、大津の方法を局所的に適用することが開発されている<ref name="jianzhuang1991automatic" />。

さらに、大津の方法の数学的根拠は、画像のヒストグラムを分散が等しくサイズが等しい2つの[[正規分布]]の混合としてモデル化する<ref name=":0">{{Cite journal |last=Kurita |first=T. |last2=Otsu |first2=N. |last3=Abdelmalek |first3=N. |date=October 1992 |title=Maximum likelihood thresholding based on population mixture models |url=https://doi.org/10.1016/0031-3203(92)90024-d |journal=Pattern Recognition |volume=25 |issue=10 |pages=1231-1240 |doi=10.1016/0031-3203(92)90024-d |issn=0031-3203}}</ref>。しかし、大津のしきい値処理はこれらの仮定を満たさない場合でも満足のいく結果をもたらす可能性がある。同様に[[統計検定]]（大津の方法が密接に関連する<ref>{{Cite journal |last=Jing-Hao Xue |last2=Titterington |first2=D. M. |date=August 2011 |title=$t$-Tests, $F$-Tests and Otsu's Methods for Image Thresholding |url=https://doi.org/10.1109/tip.2011.2114358 |journal=IEEE Transactions on Image Processing |volume=20 |issue=8 |pages=2392-2396 |doi=10.1109/tip.2011.2114358 |issn=1057-7149}}</ref>）は、動作の仮定が完全に満たされていない場合でも正しく実行される。しかし、これらの仮定から生じる深刻な逸脱を無くすために、Kittler-Illingworthの方法などの大津の方法のいくつかのバリエーションが提案されている<ref name=":0" /><ref name=":1">{{Cite journal |last=Kittler |first=J. |last2=Illingworth |first2=J. |date=1986-01-01 |title=Minimum error thresholding |url=https://www.sciencedirect.com/science/article/pii/0031320386900300 |journal=Pattern Recognition |language=en |volume=19 |issue=1 |pages=41-47 |doi=10.1016/0031-3203(86)90030-0 |issn=0031-3203}}</ref>。

=== ノイズの多い画像に対するバリエーション ===
大津の方法の限界に対処するためにさまざまな拡張が開発されてきた。人気のある拡張の1つは、ノイズの多い画像のオブジェクトセグメンテーションのタスクに適した'''2次元大津の方法'''である。ここでは、特定のピクセルの強度値をそのすぐ近傍の平均強度と比較することでセグメンテーション結果を改良する<ref name=jianzhuang1991automatic>{{cite journal
| author = Jianzhuang, Liu and Wenqing, Li and Yupeng, Tian
| year = 1991
| title = Automatic thresholding of gray-level pictures using two-dimension Otsu method
| pages = 325-327
| journal = Circuits and Systems, 1991. Conference Proceedings, China., 1991 International Conference on
}}</ref>。

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

:<math>P_{ij} = \frac{f_{ij}} N, \qquad \sum_{i=0}^{L-1}\sum_{j=0}^{L-1} P_{ij}=1</math>

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

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

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

:<math>
\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}
</math>

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

:<math>\omega_0+\omega_1 \cong 1</math>
:<math>\omega_0\mu_0+\omega_1\mu_1 \cong \mu_T</math>

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

:<math>S_b = \sum_{k=0}^1\omega_k[(\mu_k-\mu_T)(\mu_k-\mu_T)^T]</math>

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

:<math>
\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}
</math>

ここで

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

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

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

<syntaxhighlight lang="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;
</syntaxhighlight>

<math>\operatorname{tr}(S_b)</math>を評価するために、高速再起動的プログラミングアルゴリズムを使用して時間パフォーマンスを向上させることができることに留意してください<ref name=zhang2008image>{{cite journal|author1=Zhang, Jun |author2=Hu, Jinglu |name-list-style=amp | year = 2008 |title = Image segmentation based on 2D Otsu method with histogram analysis |volume = 6 |pages = 105-108 |journal = Computer Science and Software Engineering, 2008 International Conference on |doi=10.1109/CSSE.2008.206 |isbn=978-0-7695-3336-0 }}</ref>。しかし、動的プログラミングのアプローチを使用しても、2次元大津の方法は依然として時間計算量が非常に複雑である。したがって、計算コストを削減するために多くの研究が行われてきた<ref name=zhu2009fast>{{cite journal| author = Zhu, Ningbo and Wang, Gang and Yang, Gaobo and Dai, Weiming| year = 2009| title = A fast 2d otsu thresholding algorithm based on improved histogram| pages = 1-5| journal = Pattern Recognition, 2009. CCPR 2009. Chinese Conference on}}</ref>。

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

[[:en:Summed-area table]]参照

==== Matlabでの実装 ====
関数の入力と出力

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

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

<code>threshold</code>は取得されたしきい値である。

<syntaxhighlight lang="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
</syntaxhighlight>

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

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

==出典==
{{Reflist}}

==外部リンク==
* [https://github.com/cbhushan/script-collection/blob/master/GIMP-plugins/otsu-threshold.scm Implementation of Otsu's thresholding method] as [[GIMP]]-plugin using Script-Fu (a [[Scheme (programming language)|Scheme]]-based language)
* [http://homepages.inf.ed.ac.uk/rbf/CVonline/LOCAL_COPIES/MORSE/threshold.pdf Lecture notes on thresholding] – covers the Otsu method
* [http://rsb.info.nih.gov/ij/plugins/otsu-thresholding.html A plugin for ImageJ] using Otsu's method to do the threshold
* [http://www.labbookpages.co.uk/software/imgProc/otsuThreshold.html A full explanation of Otsu's method] with a working example and Java implementation
* [http://www.itk.org/Doxygen/html/classitk_1_1OtsuThresholdImageFilter.html Implementation of Otsu's method] in [[Insight Segmentation and Registration Toolkit|ITK]]
* [http://www.codeproject.com/KB/graphics/OtsuSharp.aspx Otsu Thresholding in C#] – a straightforward C# implementation with explanation
* [http://www.mathworks.com/discovery/image-thresholding.html Otsu's method using MATLAB]
* [https://scikit-image.org/docs/stable/api/skimage.filters.html#threshold-otsu Otsu Thresholding with scikit-image in Python]

{{DEFAULTSORT:おおつのにちかほう}}
[[Category:画像処理]]
[[Category:統計的偏差と分散]]