SIFT ~ 実装編(全体構成)

SIFTの理論と実装のエントリ,どうやってまとめようか悩んでたんですが,ちょっと絵を書き始めると時間がかかりすぎて前に進めなくなってしまうので,ひとまず大雑把にまとめておくことにしました.コードはだいぶリファクタして読みやすく(自分としては...)なったはずなので,コードを追いかけてもらうとSIFTの実装がどんな感じになっているかはざっくりと解ると思います.一通りSFMに目処がついた段階で,各エントリをリライトできればと思います.

あと,ブログの最後に参考文献のリンクを貼っておいたのですが,中部大の藤吉先生が作られているスライドがとっても分かりやすかったので,かっちりとした理論の部分はそちらを参考にしていただければと思います.

X.全体構成 特徴点抽出

OpenCVの cv::xfeatures2d::SIFT クラスをそのままパクってリファクタしたのですが,論文の構成にできるだけ近づくようにリファクタしました.論文の流れと一致するよう,関数を「特徴点検出」と「記述子作成」の2つに分けたのですが,ここでは特徴点検出の部分を見ています.

1.CreateInitialImage

処理の前段階として,
① 入力画像を Float 型 グレー画像へ変換.
② 入力サイズの2倍に画像を拡大.この画像が First Octave (first_octave_idx =-1) となります.
③ 2で拡大した画像に対して,ガウシアンフィルタを適用します.
を実施します.

2.Octave 数の計算.

つぎに,必要となる Octave 数を計算します.処理画像(①で前処理された画像)に対して,高さと幅の大きい方を 2 で割っていったときに,一辺の長さが 4pix を下回らないところまで画像を縮小していきますが,この値が Octave 数になります.この値は3で計算される Gaussian Pyramid の段数と同じになります.

3.BuildGaussianPyramid

画像ピラミッドを作成します.ここでは処理画像をピラミッドの底辺として,Octave 数段のピラミッドを作成します.

4.BuildDifferenceOfGaussianPyramid

3で生成した Gaussian Pyramid の各 Octave の画像間で差分を計算し,Difference of Gaussian Pyramid を作成します.

5.FindScaleSpaceExtrema

4で生成した DoG 画像の極値を取得しキーポイント候補点を取得します.候補点取得後,サブピクセル推定をします.

6.RemoveDuplicateSorted

重複しているキーポイントを削除します.

7.RetainBest

特徴点の検出上限数が指定されている場合,スコアが低いものを取り除いて上限数に収まるようにします.

8.CompensateInitialScalingToOriginalImage

1のステップで入力画像サイズを2倍にされているので,計算されたキーポイントのスケール,位置,サイズを求まった値に対して半分にします.

該当部コード
void SIFT::Detect(cv::Mat& image, std::vector<SiftKeyPoint>& keypoints) {
  int first_octave_idx = -1;

  // 1. Initialize image.
  bool double_image_size = true;
  Mat base_image = CreateInitialImage(image, double_image_size, (float)sigma_);

  // 2. Compute number of octave that is solely computed from image size.
  double log2_num = std::log((double)std::min(base_image.cols, base_image.rows)) / std::log(2.);
  int num_octaves = cvRound(log2_num - 2) + 1;

  // 3. Build Gaussian pyramid.
  std::vector<Mat> gaussian_pyr;
  BuildGaussianPyramid(base_image, gaussian_pyr, num_octaves);

  // 4. Build difference of Gaussian pyramid.
  std::vector<Mat> diff_of_gaussian_pyr;
  BuildDifferenceOfGaussianPyramid(gaussian_pyr, diff_of_gaussian_pyr);

  // 5. Find scale space extrema in difference of Gaussian pyramid.
  FindScaleSpaceExtrema(gaussian_pyr, diff_of_gaussian_pyr, keypoints);

  // 6. Remove duplicated key point.
  RemoveDuplicateSorted(keypoints);

  // 7. If maximum number supeciried.
  if (num_features_ > 0) {
    RetainBest(keypoints, num_features_);
  }

  // 8. Compensating doubling of iamge at the beginning.
  CompensateInitialScalingToOriginalImage(keypoints, first_octave_idx);
}

X.全体構成 記述子計算

次に,記述子計算部分の構成です.記述子は一つ一つのキーポイントに対して計算するため,for ループの入れ子構造になってますが,本質的な計算をしているのは下記の部分になります.

1.キーポイント周辺領域計算.

特徴点抽出の部分で抽出されたキーポイントに対して,記述子を計算するための周辺領域サイズを計算します.

2.必要なバッファ準備.

ヒストグラム計算に必要なバッファを確保します.必要なものとしては,「Mag : 勾配強度」「Ori : 勾配角度」「W : 重み(中心点からの距離による)」「RBin : 最終記述子の列方向成分」「CBin : 最終記述子の行方向成分」「hist : ヒストグラム計算領域」の6つです.

3.ComputeGradientForPrePatch

キーポイント周辺領域を分割した小領域に対して,勾配,勾配方向を計算します.

4.UpdateHistogramBasedOnTriLinearInterpolation

3では小領域に対する勾配,勾配方向を計算しましたが,これを記述子に統合します.記述子は16領域(=4x4)からなり,一つ一つの領域に8ビンのヒストグラムをもちます.小領域の計算結果がそれぞれこの16ヒストグラムの各ビンに振り分けられていきます.

5.RefineHistogram

生成されたヒストグラムを正規化し,結果を出力バッファにコピーします.

該当部コード
// Compute sift descriptor.
static void ComputeSIFTDescriptor(const Mat& img, Point2f ptf, float angle, float scl,
                                  int final_patch_width, int final_ori_hist_bins, float* dst) {
  // 1. Patch size calculation.
  int radius, pre_patch_size, final_hist_size;
  float pre_patch_radius;
  {
    pre_patch_radius = Const::SIFT_DESCR_SCL_FCTR * scl;
    // Compute radius.
    float sqrt2 = 1.4142135623730951f;
    radius = cvRound(pre_patch_radius * sqrt2 * (final_patch_width + 1) * 0.5f);
    // Clip the radius to the diagonal of the image to avoid autobuffer too large exception
    radius =
        std::min(radius, (int)sqrt(((double)img.cols) * img.cols + ((double)img.rows) * img.rows));
    pre_patch_size = (radius * 2 + 1) * (radius * 2 + 1);
    final_hist_size = (final_patch_width + 2) * (final_patch_width + 2) * (final_ori_hist_bins + 2);
  }

  // 2. Prepare buffer.
  std::vector<float> Mag(pre_patch_size, 0.0f), Ori(pre_patch_size, 0.0f), W(pre_patch_size, 0.0f),
      RBin(pre_patch_size, 0.0f), CBin(pre_patch_size, 0.0f), hist(final_hist_size, 0.0f);

  // 3. Computation for pre-patch element.
  int pix_cnt = ComputeGradientForPrePatch(
      angle, pre_patch_size, pre_patch_radius, final_patch_width, radius,
      Point2i(cvRound(ptf.x), cvRound(ptf.y)), img, RBin, CBin, Ori, Mag, W);

  // 4. Histogram update based on tri-linear interpolation.
  UpdateHistogramBasedOnTriLinearInterpolation(Mag, W, RBin, CBin, Ori, pix_cnt, angle,
                                               final_ori_hist_bins, final_patch_width, hist);

  // 5. Refine histogram.
  RefineHistogram(final_patch_width, final_ori_hist_bins, hist, dst);
}

X.実装

論文の構成にできるだけ合うようにリファクタリングしたので,もしSIFTの勉強を急遽する必要が出てきた方は下記のコードを追っかけると少しは参考になるかもです.一応,OpenCVの元の SIFT の出力と比較しながらリファクタしたので,リファクタ前後で結果は一致してます.ただ,現状実行しても可視化結果が全く出ないです..この辺も後でどうにか...
github.com

実行方法
git clone
sh build.sh
./build/sift_test

X.参考文献

Distinctive image features from scale-invariant keypoints

http://www.cs.unc.edu/~lazebnik/spring11/lec08_blob.pdf

中部大 藤吉先生の説明スライドですが,めっちゃわかりやすいです.

www.slideshare.net