大津の二値化による境界値の決定

マーチング・キューブ法でボリュームデータをメッシュにする際には、ボリューム内のCT値が一定の場所を結んでできる曲面を抽出する。この曲面を等値面と呼ぶ。

一つの材質から構成される物体のCT値は空気の部分と物体の部分の2つに分割でき、それらの場所のCT値はかなり離れた値をとっていることが期待される。このような2つの分布を単一の閾値により分割する際に閾値の値を分布の情報から自動決定する方法として大津の二値化が広く知られている。

大津の二値化の概要

大津の二値化では、とある閾値の候補を用いてデータを2つのグループに分けたときに、グループ内の平均分散 $\sigma_W^2$ (下付きの$W$はwithinの意味)とグループ間の分散$\sigma_B^2$ (下付きの$B$はbetweenの意味)を以下のように定義する。

\[\sigma_W^2 = \omega_0 \sigma_0^2 + \omega_1 \sigma_1^2\] \[\sigma_B^2 = \omega_0 (\mu_0 - \mu)^2 + \omega_1 (\mu_1 - \mu)^2\]

上記の式において、$\mu$はデータ全体の平均、$\mu_0, \mu_1$は各グループでの平均、$\sigma_0^2, \sigma_1^2$は各グループの分散、$\omega_0, \omega_1$は各グループに含まれるデータの割合を表す。それぞれは値$i$を取る確率を$p(i)$, データ長さを$N$, 閾値を$T$として以下のように定義される。

\[\mu = \sum_{i=0}^{N-1} i p(i), \quad \mu_0 = \frac{\displaystyle \sum_{i=0}^{T-1} i p(i)}{\displaystyle \sum_{i=0}^{T-1} p(i)}, \quad \mu_1 = \frac{\displaystyle \sum_{i=T}^{N-1} i p(i)}{\displaystyle \sum_{i=T}^{N-1} p(i)}\] \[\sigma_0^2 = \frac{\displaystyle \sum_{i=0}^{T-1} (\mu_0 - i)^2 p(i)}{\displaystyle \sum_{i=0}^{T-1} p(i)}, \quad \sigma_1^2 = \frac{\displaystyle \sum_{i=T}^{N-1} (\mu_1 - i)^2 p(i)}{\displaystyle \sum_{i=T}^{N-1} p(i)}\] \[\omega_0 = \sum_{i=0}^{T-1} p(i), \quad \omega_1 = \sum_{i=T}^{N-1} p(i)\]

これらを用いると、2つのグループの分散度$\lambda$が、グループ内分散とグループ間分散の比として以下のように表せる。

\[\lambda = \frac{\sigma_B^2}{\sigma_W^2}\]

ただし、実際には$\sigma_B^2 = \sigma^2 - \sigma_W^2$であることから、$\sigma_B^2$が最大となるような閾値を求めれば十分で、さらにグループ間分散は平均だけを使って以下のように書き直せる。

\[\sigma_B^2 = \omega_0 \omega_1 (\mu_0 - \mu_1)^2\]

この値をすべての閾値の候補に対して計算し、最も大きな値を取る(=データが最もはっきり分割できる)閾値を求めるのが大津の二値化のアルゴリズムである。

大津の二値化の実装

上記の式に従い、平均と分散を求めようとすると、プログラムとしての効率はどうだろうか？

今、ボリュームデータに含まれているデータはそのCT値の大きさを元にヒストグラム化されているとしよう。データの範囲が$[0, N-1]$とすれば、ヒストグラムとは各 $i \in [0, N-1]$につき、データの出現確率$p(i)$が求まっていることに相当する。

このヒストグラムの値を用いると値の平均ならびに分散は以下のように計算できる。

\[\mu = \sum_{i=0}^{N-1} i p(i)\]

同様に、分散は以下のように計算できる。

\[\sigma^2 = \sum_{i=0}^{N-1} (i - \mu)^2 p(i)\]

これはグループ内の平均や分散についても、同様なので、各グループで平均や分散を求める操作の計算量は$O(N)$と見積もられる。これは上記の式で定義される$\mu_0, \mu_1$を求める場合も同様で、オーダーとしての計算量は$O(N)$である。

今、閾値の候補は$N$個あるから、このままだと計算量は$O(N^2)$になりそうだ。しかし、CT値が16bitの符号なし整数、すなわち0-65535の値で表されていることを考えると$N^2$は$10^9$程度のオーダーとなってしまい、計算量が大きすぎる。

余談になるが、現代のCPUは3GHz程度のクロック数なので、単純な計算が1秒間に$3 \times 10^9$できるが、実際の計算では目安として計算量が$10^7$のオーダーを超えると一瞬では計算が終わらなくなる。

効率的な実装

上記の通り、大津の二値化ではグループ間分散が最大となる閾値を選べばよく、この値は各グループ内平均 $\mu_0, \mu_1$から計算できる。繰り返しとなるが、閾値$T$については、これらのグループ内平均は次のように表せる。

\[\begin{aligned} \mu_0 &= \frac{1}{\omega_0} \sum_{i=0}^{T-1} i p(i) \\ \mu_1 &= \frac{1}{\omega_1} \sum_{i=T}^{N-1} i p(i), \\ \end{aligned}\]

この式をよく見ると、Tが1増えたときには、以下の操作を行うことで、平均の値が簡単に更新できることが分かる。

• $\omega_0$を$p(T)$増やし、$\omega_1$を$p(T)$減らす。
• $\Sigma$の中に含まれる $i p(i)$を$\mu_1$の側から$\mu_0$の側に移項する。

この操作をすれば、いちいち、各閾値$T$に対してヒストグラムの値を足し上げる操作をする必要がないので、非常に効率的であることが分かる。上記の2つの操作はデータの値の範囲に寄らず、同様の操作は各グループの分散を求めるときにも使えるので、結果として大津の二値化の計算量が$O(N)$まで減少する。

ここまでの議論をまとめると、効率的な大津の二値化の実装では、はじめに$T=0$のときのグループ1に対する平均と出現頻度$\omega_1 = 1$を求めておき、閾値を増やすごとに、上記の移項の操作を行いながら、最大の分散度を与える閾値を探索すれば良い。

サンプルプログラムのmcubes.cpp内にある関数 getThresholdOtsu に大津の二値化を実装してみよう。