MENU

C++でOpenCVを使って判別分析法(大津の2値化)を自作してみる

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

判別分析法とは

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

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

opencv-discriminant-analysis-method-otsu-cpp-1

閾値 $ t $ で2値化するとき、クラス1の画素数を $ w_1 $ 、平均(輝度値の合計÷画素数)を $ m_1 $ 、分散を $ \sigma_1^2 $ 、同様にしてクラス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);

こんな感じで書けば判別分析法で2値化できます。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;  // クラス1の画素数
    int w2 = 0;  // クラス2の画素数
    long sum1 = 0;  // クラス1の平均を出すための合計値
    long sum2 = 0;  // クラス2の平均を出すための合計値
    double m1 = 0.0;  // クラス1の平均
    double m2 = 0.0;  // クラス2の平均
    
    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の値を使って2値化 */
  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 関数内で、両者の戻り値を表示し一致していれば自作関数が正しく動作していることになります。