再び、判別分析法による閾値の自動計算について(2)

2012 / 01 / 20 by
Filed under: AS 
Bookmark this on Delicious
[`livedoor` not found]
[`yahoo` not found]

最近ネットで見つけた「大津の二値化画像」というページに書いてある判別分析法による閾値の自動計算方法が、著しく少ないループ回数で済んでいる、という話の続き。前回はこちら

ステップ 3. のループでは2つの新たな配列を生成しています。w_mo は画素数(厳密に言うと画素数の割合)を、u_mo は濃度を、ヒストグラムから計算して格納していますが、単純に各諧調の値を格納しているわけではありません。1つ前の要素と当該諧調の値を足合わせたものを格納しています。

ということはΣi=0Nの演算をおこなっているわけですね(N は配列のインデックス番号)。

例えばインデックス 10 の場所には、Σi=010、つまりインデックス 0 ~10 までのすべての値を加算した値が格納され、インデックス 20 の場所には、Σi=020、つまりインデックス 0 ~20 までのすべての値を加算した値が格納され…… といった具合。

ここで再び、以前「二値化(3) ~判別分析法による閾値の自動計算~」で組んだ、とんでもないループ回数のコードを見てみます。

二重ループでは外側のインデックスの値を境にして、内側のループを2分割してます。
例えば外側のインデックスが55の場合、内側のループでは256回のループを 0~54 と 55~255 の2つに分割、外側のインデックスが127の場合、内側のループでは256回のループを 0~126 と 127~255 の2つに分割、といった具合。

で、内側の2分割されたそれぞれのループの中では、Σi=054とΣi=55255やΣi=0126とΣi=127255などの計算をおこなっているわけです。

つまり「大津の二値化画像」のステップ 3. では、二重ループの内側のループで内側で2分割されるもののうち、先のループで計算される値を、あらかじめ別の配列として用意しているわけです。

で、内側の2分割されたループの後のループで求める値は、w_mo[255]、u_mo[255] から、2分割されたループのうち、前のループで得られた値を引けば得られる、と。

なるほど! これなら外側のループが回るたびに毎回毎回おこなっていた内側のループを、別途1回だけで済ませられるわけだ。あったまいい。

でもステップ 4. の class_bunsan を求める式の原理がよく分らん。
いろいろ調べたけど、こういう風に濃度と画素数から分散を求める公式を使っている事例が他に見つかりませんでした。

ところで以前「二値化(3) ~判別分析法による閾値の自動計算~」でコードを組んだ時、参考にした書籍がもう一冊ありました。
C言語で学ぶ実践画像処理」という本ですが、それを読み返したらば、何となく似たような計算をしていることに気づきました。P25~26 の説明によると、クラス内・クラス間という二つの分散ではなく、以下の一つの分散を計算すれば良いとあります。


ρ^2 = (n1 * (m1 - m0)^2 + n2 * (m2 - m0)^2) / (n1 + n2)
 ρ : 分散
 m0, m1, m2 : 全体の濃度、物体1の濃度、物体2の濃度
 n1, n2 : 物体1の画素数、物体2の画素数

この分散を256階調すべてについて求め、最大になるときの階調を閾値とする、とのこと。
以前のやり方だと m1, m2, n1, n2 を256の各諧調ごとにおこなうので 256 * 256 = 65,536 回のループが必要になります。
が、「大津の二値化画像」のクレバーなステップ 3. を使えば 256 * 256 = 65,536 回ではなく 256 + 256 = 512 回にまでループ回数を圧縮できます。

何と 130,048 回から 512 回に減らすことができました。1 / 254 という驚異の数値!

C言語で学ぶ実践画像処理」と「大津の二値化画像」を参考に書いた判別分析法による閾値の自動計算するクラスを以下に示します。

package {
  import flash.display.BitmapData;
  /**
   * 判別分析法による閾値の計算
   */
  public class DiscriminantAnalysis {
    /**
     * 判別分析法による閾値の計算
     * @param  bmd  対象イメージ(グレイスケール化されていること)
     * @return  計算によって求められた閾値(最大分散)
     */
    static public function calculate(bmd:BitmapData):uint {
      // n : 画素数、 m : 濃度、 p : 分散
      // S0 : 全体、 S1 : 指定閾値より小さい値の領域、 S2 : 指定閾値より大きい値の領域
      
      // ヒストグラム取得(RGBA のうち A 以外で便宜上 R を使用)
      var hist:Vector.<Number> = bmd.histogram()[0];
      
      // 各階調における画素数と濃度の総和を求める
      var sumN:Vector.<uint> = new Vector.<uint>(256, true);
        // あるインデックスに格納された値はそれ以下の全階調の画素数をすべて足したもの
        // 例えば sumN[5] に格納された値は hist[0] ~ hist[5] までのすべての値を足したもの
      var sumM:Vector.<uint> = new Vector.<uint>(256, true);
        // sumN と同様にデータを格納するが、その値はヒストグラムから求められる濃度(画素数 * 階調)
        // 例えば sumM[5] に格納された値は hist[0] * 0 ~ hist[5] * 5 までのすべての値を足したもの
      var i:int = 0;
      sumN[i] = hist[i];
      sumM[i] = hist[i] * i;
      for (i = 1; i < 256; i++) {
        sumN[i] = sumN[i - 1] + hist[i];
        sumM[i] = sumM[i - 1] + hist[i] * i;
      }
      
      // 全体の各値を求める
      var n0:uint = sumN[255];  // 画素数(= bmd.width * bmd.height)
      var m0:uint = sumM[255];  // 濃度
      
      // 各階調の値を閾値としたときの分散の計算
      var maxP:Number    = 0;  // 最大分散
      var threshold:uint = 0;    // 求める閾値
      for (i = 0; i < 256; i++) {
        // 指定閾値より小さい値の領域
        var n1:uint   = sumN[i];  // 画素数
        var m1:Number = sumM[i];  // 濃度
        // 指定閾値より大きい値の領域
        var n2:uint   = n0 - n1;  // 画素数
        var m2:Number = sumM[255] - m1;  // 濃度
        // 領域の濃度の平均
        if (n1 == 0) { m1 = 0; } else { m1 /= n1; }  // 指定閾値より小さい値の領域
        if (n2 == 0) { m2 = 0; } else { m2 /= n2; }  // 指定閾値より大きい値の領域
        
        // 分散計算
        var val1:Number = m1 - m0;
        var val2:Number = m2 - m0;
        var p:Number = (n1 * val1 * val1 + n2 * val2 * val2) / n0;
        if (p > maxP) {
          maxP      = p;
          threshold = i;
        }
      }
      return threshold;
    }
  }
}


Comments

Tell me what you're thinking...
and oh, if you want a pic to show with your comment, go get a gravatar!





WP-SpamFree by Pole Position Marketing