「大津の手法」を実装しました
どうも, Boltzmanです. 魚眼らしいことがまだ全然できていないAR, 3回目です.
今回の目標は「大津の手法」(判別分析法)による閾値自動決定の実装です. 大津の手法を紹介しているサイトはたくさんありますが, ここでも一通り紹介したいと思います. 正直書くことがない.
大津の手法は「分離度」と呼ばれる値が最大になるような閾値を採用する方法で, 分離度は[\frac{\sigma_b^2}{\sigma_w^2}]で表されます. ここで(\sigma_b)はクラス間分散, (\sigma_w)はクラス内分散です.
分散は要素がどれくらい散らばっているかの指標です.
上の図では, 右側のほうが左側に比べて要素が散らばっているため, 右のほうが分散の値は大きくなります. 分散(\sigma)は, 要素(\boldsymbol{x}i)とそれらの平均(\boldsymbol{\mu})を用いて\begin{equation}\sigma^2=\sum{i}(\boldsymbol{x}_i-\boldsymbol{\mu})^2\end{equation}と表されます. つまり, 分散は要素の平均と各要素との距離の和と考えることができます.
ここで要素を複数のクラスに分けます.
上のデータを見るとなんとなく二つに分けられそうな気がします.
こんな風に. このようにして括った二つのクラスについて, 分散を考えます.
まずクラス内分散は各クラス内の要素の分散です. つまり, \begin{equation}\sigma_w^2=\sum_{k}\sum_{\boldsymbol{x}_i\in\chi_k}\frac{\omega_k}{\omega}(\boldsymbol{x}_i-\boldsymbol{\mu}_k)^2\end{equation}こんな感じですかね. (\boldsymbol{\mu}_k)はクラス内での平均, (\omega_k)はクラス内の要素数, (\omega)は全体の要素数を表します.
クラス間分散は, 各クラスがどれくらい散らばっているかを示します. つまり, \begin{equation}\sigma_b=\sum_{k}\frac{\omega_k}{\omega}(\boldsymbol{\mu}-\boldsymbol{\mu}_k)^2\end{equation}こうですかね. (\boldsymbol{\mu})は全体の平均を表します.
なんか, ちゃんと統計やらなきゃだめだなって思います. 今回は実装だけなので原理はかなりスルーしてます. 適当なこと書いてたらすいません.
で, 分離度はクラス内分散(\omega_w)が小さいほど, クラス間分散(\omega_b)が大きいほど大きくなります. つまり, 各クラスを見たときに要素が集まっているほど, クラス同士を見たときに各クラスが散らばっているほど大きくなります.
ここで, クラス内分散(\sigma_w)とクラス間分散(\sigma_b), それに全分散(\sigma)について, 以下の関係式を導入します. [\sigma^2=\sigma_w^2+\sigma_b^2]この証明ができなくて困っています. どのサイト見ても自明っぽく書いてあるので, たぶん私が何か重要な公式を知らないんじゃないだろうか. この関係式を認めるとすると, 分離度は以下のように書きなおせます. [\frac{\sigma_b^2}{\sigma^2-\sigma_b^2}]よって, クラス間分散が最大になる時, 分離度は最大になります. なぜクラス間分散のほうを残したのかというと, クラス内分散のほうは計算の時に各クラスでの分散を求める必要があって面倒だからです.
今回はこれを使って, 濃淡のヒストグラムを二つに分割します.
閾値でヒストグラムを二つのクラスに分けた時の分離度を考えます. クラス間分散は上の定義から, \begin{equation}\sigma_b^2=\frac{\omega_1(\mu_1-\mu)^2+\omega_2(\mu_2-\mu)^2}{\omega}\end{equation}ここで, [\mu=\frac{\omega_1\mu_1+\omega_2\mu_2}{\omega}]より, \begin{eqnarray}\mu_1-\mu &=&\mu_1-\frac{\omega_1\mu_1+\omega_2\mu_2}{\omega}\&=&\frac{\omega_1(\mu_1-\mu_2)}{\omega}\end{eqnarray}となるので, \begin{equation}\sigma_b^2=\frac{\omega_1\omega_2(\mu_1-\mu_2)^2}{\omega^2}\end{equation}と書き直せます. (\omega)はクラスの分け方によらず一定なので, (\omega_1\omega_2(\mu_1-\mu_2)^2)が最大となる時(\sigma_b)も最大となり, 分離度が最大になります.
上の図で, ヒストグラムの式を[y=f(x)(y\in\mathbb{Z}:要素数, x\in\mathbb{Z}:グレースケールの値(0\leq x\leq 255))]とすると, 平均は\begin{equation}\mu_k=\frac{1}{\omega_k}\sum_{x_i\in \chi_k}x_if(x_i)\end{equation}であたえられます. (平均の定義を勘違いしてたのは私だけでいい)
結局コードは以下のようになりました.
</p>
<p>public int getThresh()<br />
{<br />
int th = 0; //ループ用の閾値<br />
double num1, num2; //画素数<br />
double sum1, sum2; //和<br />
double m1, m2; //平均<br />
double va_b; //クラス間分散の二乗<br />
double max = 0; //va_bの最大値を保存しておく変数</p>
<p> //ヒストグラムを生成していなければ弾く<br />
if (hist == null) return -1;</p>
<p> //(th=0)の処理<br />
num1 = hist[0];<br />
sum1 = hist[0] * 0;<br />
m1 = num1 != 0 ? sum1 / num1 : 0;</p>
<p> num2 = 0;<br />
sum2 = 0;<br />
m2 = 0;<br />
for (int i = 1; i &lt; 256; i++)<br />
{<br />
num2 += hist[i];<br />
sum2 += hist[i] * i;<br />
}<br />
m2 = num2 != 0 ? sum2 / num2 : 0;</p>
<p> va_b = (num1 * num2 * (m1 – m2) * (m1 – m2));<br />
max = va_b;<br />
//Debug.Log("0:"+va_b);</p>
<p> //(th=1)以降の処理<br />
for (th = 1; th &lt; 255; th++)<br />
{<br />
num1 += hist[th];<br />
num2 -= hist[th];<br />
sum1 += hist[th] * th;<br />
sum2 -= hist[th] * th;<br />
m1 = num1 != 0 ? sum1 / num1 : 0;<br />
m2 = num2 != 0 ? sum2 / num2 : 0;<br />
va_b = (num1 * num2 * (m1 – m2) * (m1 – m2));</p>
<p> //Debug.Log(th+":"+va_b);<br />
if (max < va_b)<br />
{<br />
thresh = th;<br />
max = va_b;<br />
}<br />
}</p>
<p> return thresh;<br />
}</p>
<p>
これで一応動いていますが, 原理が完全には腑に落ちていないせいか, どうにも達成感がありません.
今回はここまでです. 次回はシェーダーを勉強しないといけないと思うので, 遅れると思います. では.