C++でOpenCVを使って判別分析法（大津の２値化）を自作してみる

OpenCVの関数に頼っているとアルゴリズムなんて分からなくても使えてしまうので、画像処理の勉強をしている以上、OpenCVの関数を使うときはその仕組みを理解してから使おうと思います。今回は、２値化アルゴリズムの紹介です。

判別分析法とは

大津の２値化とも呼ばれ、画像を２値化する手段として最も一般的に利用されているアルゴリズムです。他のアルゴリズムと違い、閾値を自動的に算出するため、非常に便利な２値化法です。

クラス間分散とクラス内分散により求められる分離度が最大値となる閾値で２値化します。

閾値 $ t $ で２値化するとき、クラス１の画素数を $ w_1 $ 、平均（輝度値の合計÷画素数）を $ m_1 $ 、分散を $ \sigma_1^2 $ 、同様にしてクラス２の画素数、平均、分散をそれぞれ、$ w_2 $ 、$ m_2 $ 、$\sigma_2^2 $ とします。また、画像全体の画素数を $ w_t $ 、平均を $ m_t $ 、分散を $ \sigma_t^2 $ としたとき、クラス内分散 $ \sigma_w^2 $ は、
$$ \sigma_w^2 = \frac{w_1 \sigma_1^2 + w_2 \sigma_2^2}{w_1 + w_2} $$

クラス間分散 $ \sigma_b^2 $ は、
$$ \begin{eqnarray*}
\sigma_b^2 &=& \frac{w_1 (m_1 – m_t)^2 + w_2 (m_2 – m_t)^2}{w_1 + w_2} \\
&=& \frac{w_1 w_2 (m_1 – m_2)^2}{(w_1 + w_2)^2}
\end{eqnarray*} $$

ここで、$ \sigma_t^2 = \sigma_w^2 + \sigma_b^2 $ であるから、分離度（判別比）は、
$$ \frac{\sigma_b^2}{\sigma_w^2} = \frac{\sigma_b^2}{\sigma_t^2 – \sigma_b^2} $$

また、$ \sigma_t^2 $ は常に一定なので $ \sigma_b^2 $ が最大となる閾値を求めればよいとわかる。さらに、$ \sigma_b^2 $ の式に着目すると、分母の $ (w_1 + w_2)^2 $ も常に一定なので、分子の $ w_1 w_2 (m_1 – m_2)^2 $ が最大となるときの閾値 $ t $ を求めればよい。

ソースコード

OpenCVの関数を使うと

cv::threshold(入力画像, 出力画像, 0, 255, cv::THRESH_BINARY|cv::THRESH_OTSU);

こんな感じで書けば判別分析法で２値化できます。cv::THRESH_OTSU というところで大津法を指定していますね。

自作

#include <iostream>
#include <vector>
#include <opencv2/imgproc.hpp>
#include <opencv2/highgui.hpp>

int discriminantAnalysis(cv::Mat_<uchar> src, cv::Mat_<uchar>& dst){
  /* ヒストグラム作成 */
  std::vector<int> hist(256, 0);  // 0-255の256段階のヒストグラム（要素数256、全て0で初期化）
  for (int y = 0; y < src.rows; ++y){
    for (int x = 0; x < src.cols; ++x){
      hist[static_cast<int>(src(y, x))]++;  // 輝度値を集計
    }
  }
  
  /* 判別分析法 */
  int t = 0;  // 閾値
  double max = 0.0;  // w1 * w2 * (m1 - m2)^2 の最大値
  
  for (int i = 0; i < 256; ++i){
    int w1 = 0;  // クラス１の画素数
    int w2 = 0;  // クラス２の画素数
    long sum1 = 0;  // クラス１の平均を出すための合計値
    long sum2 = 0;  // クラス２の平均を出すための合計値
    double m1 = 0.0;  // クラス１の平均
    double m2 = 0.0;  // クラス２の平均
    
    for (int j = 0; j <= i; ++j){
      w1 += hist[j];
      sum1 += j * hist[j];
    }
    
    for (int j = i+1; j < 256; ++j){
      w2 += hist[j];
      sum2 += j * hist[j];
    }
    
    if (w1)
      m1 = (double)sum1 / w1;
    
    if (w2)
      m2 = (double)sum2 / w2;
    
    double tmp = ((double)w1 * w2 * (m1 - m2) * (m1 - m2));
    
    if (tmp > max){
      max = tmp;
      t = i;
    }
  }
  
  /* tの値を使って２値化 */
  for (int y = 0; y < src.rows; ++y){
    for (int x = 0; x < src.cols; ++x){
      if (src(y, x) < t)
        dst(y, x) = 0;
      else
        dst(y, x) = 255;
    }
  }
  
  return t;
}

int main(){
  cv::Mat_<uchar> src = cv::imread("image.jpg", 0);  // グレースケールで画像読み込み
  cv::Mat_<uchar> src2 = cv::imread("image.jpg", 0);  // グレースケールで画像読み込み
  
  int a = cv::threshold(src, src, 0, 255, cv::THRESH_BINARY|cv::THRESH_OTSU);
  int b = discriminantAnalysis(src2, src2);
  
  std::cout << a << std::endl;
  std::cout << b << std::endl;
}

cv::threshold() は戻り値として閾値 $ t $ を返すので、自作関数も閾値を返すようにしました。main 関数内で、両者の戻り値を表示し一致していれば自作関数が正しく動作していることになります。