大津の二値化法

大津の二値化法（おおつのにちかほう、Otsu's method、大津の方法などとも呼ばれる）とは、コンピュータビジョンや画像処理において、自動画像しきい値処理を行うために使用される手法。その名はにちなむ。最も単純な形式では、アルゴリズムはピクセルをforegroundとbackgroundの2つのクラスに分ける1つの強度しきい値を返す。このしきい値はクラス内の強度分散を最小化することにより、または同等にクラス間の分散を最大化することにより決定される。
大津の二値化はフィッシャーの判別分析の1次元の離散に相当するものであり、に関連しており、強度ヒストグラムで行われる大域最適なK平均法と同等である。マルチレベルのしきい値処理へ拡大することは最初の論文で説明されており、以降計算効率の高い実装が提案されている。
大津の二値化.
アルゴリズムは、2つのクラスの分散の加重和として定義されるクラス内分散を最小化するしきい値を探す。
加重formula_2とformula_3はしきい値formula_4により分けられた2つのクラスの確率であり、formula_5とformula_6はこれら2つのクラスの分散である。
クラスの確率formula_7はヒストグラムのformula_8個のビンから計算される。
2つのクラスでは、クラス内分散を最小化することはクラス間分散を最大化することと同じである。
これはクラス確率formula_11とクラス平均formula_12で表され、クラス平均formula_13、formula_14およびformula_15は、
である。以上の関係は以下のように簡単に確かめることができる。
クラス確率とクラス平均は、繰り返し計算することができる。このアイデアは効果的なアルゴリズムを生む。
MATLABまたはOctaveでの実装.
codice_1はさまざまなグレーレベルのグレースケール画像（通常は8ビット画像）の256画素のヒストグラムである。codice_2は画像のしきい値 (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に組み込み関数codice_3とcodice_4があり、それぞれ大津の手法とマルチ大津の手法(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.
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次元大津の方法である。ここでは、特定のピクセルの強度値をそのすぐ近傍の平均強度と比較することでセグメンテーション結果を改良する。
各ピクセルで、近傍の平均グレーレベル値が計算される。与えられたピクセルのグレーレベルをformula_8個の離散値に分割し、平均グレーレベルも同じformula_8個の値に分割する。その後、ピクセルのグレーレベルと近傍の平均のペアformula_27が形成される。各ペアはformula_28個の考えうる2次元ビンの1つに属する。ペアformula_27の出現の総数（頻度）formula_30を画像のピクセルの総数formula_31で割ったものは、2次元ヒストグラムで同時確率密度関数を定義する。
そして、2次元大津の方法は、2次元ヒストグラムに基づいて次のように開発されている。
2つのクラスの確率は、次のように表せる。
2つのクラスの強度平均ベクトルと合計平均ベクトルは次のように表せる。
ほとんどの場合、非対角の確率はわずかであるため、次のことを簡単に確認することができる。
クラス間離散行列は次のように定義される。
離散行列のトレースは次のように表される。
ここで
1次元の大津の方法と同様に、最適なしきい値formula_41は、formula_42を最大化することで取得される。
アルゴリズム.
formula_43とformula_4は、1次元の大津の方法と同様に繰り返し取得される。formula_43とformula_4の値は、formula_42の最大値を取得するまで変更される。つまり、
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;
formula_42を評価するために、高速再起動的プログラミングアルゴリズムを使用して時間パフォーマンスを向上させることができることに留意してください。しかし、動的プログラミングのアプローチを使用しても、2次元大津の方法は依然として時間計算量が非常に複雑である。したがって、計算コストを削減するために多くの研究が行われてきた。
合計領域テーブルを使用して3つのテーブルを作製する場合は、formula_49を合計し、formula_50を合計し、formula_51を合計し、その時のランタイムの複雑性は(O(N_pixels), O(N_bins*N_bins))の最大値である。しきい値に関して粗い解像度のみが必要な場合はN_binsを減らすことができることに注意してください。
参照
Matlabでの実装.
関数の入力と出力
codice_5は、グレースケール値と近傍平均グレースケール値のペアのformula_52の2次元ヒストグラムである。
codice_6は、所与の画像のペアの数である。これは各方向の2次元ヒストグラムのビンの数により決まる。
codice_7は取得されたしきい値である。
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のアルゴリズムよりも大津の方法を優先する方が好ましい場合がある。