推定器

■ Condensation

サンプル

Condensation

Condensationアルゴリズムを用いて画像中の赤い領域をトラッキングする

サンプルコード

#include <cv.h>
#include <highgui.h>
#include <ctype.h>
#include <math.h>

// (1)尤度の計算を行なう関数
float
calc_likelihood (IplImage * img, int x, int y)
{
  float b, g, r;
  float dist = 0.0, sigma = 50.0;

  b = img->imageData[img->widthStep * y + x * 3];	//B
  g = img->imageData[img->widthStep * y + x * 3 + 1];	//G
  r = img->imageData[img->widthStep * y + x * 3 + 2];	//R
  dist = sqrt (b * b + g * g + (255.0 - r) * (255.0 - r));

  return 1.0 / (sqrt (2.0 * CV_PI) * sigma) * expf (-dist * dist / (2.0 * sigma * sigma));
}


int
main (int argc, char **argv)
{
  int i, c;
  double w = 0.0, h = 0.0;
  CvCapture *capture = 0;
  IplImage *frame = 0;

  int n_stat = 4;
  int n_particle = 4000;
  CvConDensation *cond = 0;
  CvMat *lowerBound = 0;
  CvMat *upperBound = 0;

  int xx, yy;

  // (2)コマンド引数によって指定された番号のカメラに対するキャプチャ構造体を作成する
  if (argc == 1 || (argc == 2 && strlen (argv[1]) == 1 && isdigit (argv[1][0])))
    capture = cvCreateCameraCapture (argc == 2 ? argv[1][0] - '0' : 0);

  // (3)１フレームキャプチャし，キャプチャサイズを取得する．
  frame = cvQueryFrame (capture);
  w = frame->width;
  h = frame->height;

  cvNamedWindow ("Condensation", CV_WINDOW_AUTOSIZE);

  // (4)Condensation構造体を作成する．
  cond = cvCreateConDensation (n_stat, 0, n_particle);

  // (5)状態ベクトル各次元の取りうる最小値・最大値を指定する．
  lowerBound = cvCreateMat (4, 1, CV_32FC1);
  upperBound = cvCreateMat (4, 1, CV_32FC1);

  cvmSet (lowerBound, 0, 0, 0.0);
  cvmSet (lowerBound, 1, 0, 0.0);
  cvmSet (lowerBound, 2, 0, -10.0);
  cvmSet (lowerBound, 3, 0, -10.0);
  cvmSet (upperBound, 0, 0, w);
  cvmSet (upperBound, 1, 0, h);
  cvmSet (upperBound, 2, 0, 10.0);
  cvmSet (upperBound, 3, 0, 10.0);

  // (6)Condensation構造体を初期化する
  cvConDensInitSampleSet (cond, lowerBound, upperBound);

  // (7)ConDensationアルゴリズムにおける状態ベクトルのダイナミクスを指定する
  cond->DynamMatr[0] = 1.0;
  cond->DynamMatr[1] = 0.0;
  cond->DynamMatr[2] = 1.0;
  cond->DynamMatr[3] = 0.0;
  cond->DynamMatr[4] = 0.0;
  cond->DynamMatr[5] = 1.0;
  cond->DynamMatr[6] = 0.0;
  cond->DynamMatr[7] = 1.0;
  cond->DynamMatr[8] = 0.0;
  cond->DynamMatr[9] = 0.0;
  cond->DynamMatr[10] = 1.0;
  cond->DynamMatr[11] = 0.0;
  cond->DynamMatr[12] = 0.0;
  cond->DynamMatr[13] = 0.0;
  cond->DynamMatr[14] = 0.0;
  cond->DynamMatr[15] = 1.0;

  // (8)ノイズパラメータを再設定する．
  cvRandInit (&(cond->RandS[0]), -25, 25, (int) cvGetTickCount ());
  cvRandInit (&(cond->RandS[1]), -25, 25, (int) cvGetTickCount ());
  cvRandInit (&(cond->RandS[2]), -5, 5, (int) cvGetTickCount ());
  cvRandInit (&(cond->RandS[3]), -5, 5, (int) cvGetTickCount ());

  while (1) {
    frame = cvQueryFrame (capture);

    // (9)各パーティクルについて尤度を計算する．
    for (i = 0; i < n_particle; i++) {
      xx = (int) (cond->flSamples[i][0]);
      yy = (int) (cond->flSamples[i][1]);
      if (xx < 0 || xx >= w || yy < 0 || yy >= h) {
	cond->flConfidence[i] = 0.0;
      }
      else {
	cond->flConfidence[i] = calc_likelihood (frame, xx, yy);
	cvCircle (frame, cvPoint (xx, yy), 2, CV_RGB (0, 0, 255), -1);
      }
    }

    cvShowImage ("Condensation", frame);
    c = cvWaitKey (10);
    if (c == '\x1b')
      break;

    // (10)次のモデルの状態を推定する
    cvConDensUpdateByTime (cond);

  }

  cvDestroyWindow ("Condensation");
  cvReleaseImage (&frame);
  cvReleaseCapture (&capture);
  cvReleaseConDensation (&cond);
  cvReleaseMat (&lowerBound);
  cvReleaseMat (&upperBound);

  return 0;
}

// (1)尤度の計算を行なう関数
Condensationアルゴリズムでは，各パーティクルの尤度を与える必要がある． サンプルコードでは，各パーティクルの位置の情報から，そのピクセル値の赤らしさを 尤度として用いる事にした．直感的な分かりやすさのためにRGB表色系のまま計算を 行なう．具体的には以下のように，R(255,0,0)からのユークリッド距離 dに対して，0を平均， sigmaを分散\sigmaとして持つような正規分布を 尤度関数L(d)として定義している．

// (2)コマンド引数によって指定された番号のカメラに対する キャプチャ構造体を作成する
関数 cvCreateCameraCapture()を用いて， コマンド引数として与えられた番号のカメラに対するキャプチャ構造体 CvCaptureを作成する． 引数が指定されない場合は，デフォルトの値である"0"が利用される．

// (3)1フレームキャプチャし，キャプチャサイズ を取得する．
関数 cvQueryFrame() を呼び出し，実際に１フレー ムキャプチャを行い，取得した画像から，画像の幅wと高さhを取得する．

// (4)Condensation構造体を作成する．
関数 cvCreateConDensation() を，状態ベクトルの次元数n_statとパーティクルの数 n_particleを指定して呼び出し，Condensatopn構造体を作成する． 第２引数は，観測ベクトル数を指定する事になっているが，OpenCVの現バージョンでは 観測ベクトルは使用していないため，値には0を指定したのでかまわない．
サンプルコードでは，ダイナミクスとして以下のように等速直線運動を仮定しているので， 状態ベクトルはその位置(x,y)と各軸方向の速度(u,v)の4次元となる．

// (5)状態ベクトル各次元の取りうる最小値・最大値を指定する．
ベクトル lowerBound, upperBound をそれぞれ関数cvCreateMat()で領域確保し， 状態ベクトルの各次元の変数が取りうる値の下限値（最小値）と上限値（最大値）を設定する． この値は，初期値設定の際にはパーティクルを一様分布させるために使用される．
サンプルコードでは，位置に対して最小値0，最大値は画像サイズ(w,h)を指定，また速度に 対しては -10〜10[pixel]を指定している．

// (6)Condensation構造体を初期化する
関数 cvConDensInitSampleSet() で，構造体の初期化，初期パーティクルの状態の設定を 行なう．

// (7)ConDensationアルゴリズムにおける状態ベクトルのダイナミクスを指定する
CvConDensation構造体のメンバDynamMatr[]にダイナミクスを表現する行列要素を順に指定する．
サンプルコードでは，(4)で示した等速直線運動を以下のように行列表現した際の4X4の配列を与えている．

// (8)ノイズパラメータを再設定する．
初期化の際に，状態ベクトルの各次元の変数に対するノイズパラメータは自動的に設定される （上限値，下限値それぞれの範囲の1/5の値になる）ので，それが不都合な場合はここで再設定する．
サンプルコードでは，位置に関しては-25〜25[pixel]，速度に関しては-5〜5[pixel]の範囲のノイズ を指定している．

// (9)各パーティクルについて尤度を計算する．
ループブロック中で関数 cvQueryFrame() を呼び出し，画像を取得する． その後，各パーティクルの位置（サンプルコードでは，CvConDensation構造体のflSamples[i][0]と[i][1]） を用いてその尤度を(1)で定義した尤度関数を用いて計算する．
パーティクルの位置が画像範囲内であれば，青い丸をキャプチャした画像上に追加し，画像を描画する． "Esc"キーが押された場合は，トラッキングを終了する．

// (10)次のモデルの状態を推定する
関数cvConDensUpdateByTime()を用いて，次の時刻のモデルの状態を推定する． condensationアルゴリズムにおける，リサンプル・移動・拡散のフェーズを行なっている と考えてよい．

実行結果例