2021年と鳩の糞

あけましておめでとうございます!

年も明けた1月1日,クリスティアーノ・ロナウドばりの腹筋を維持するために,ランニングと筋トレに嫌々(笑)出かけたのですが,いつものランニングコースを走っているとなにやら上から落下物が....「まさか!」と思ったのですが,鳩の糞でした...

うーん.これはどう解釈すればいいの?と思ってたんですが,ネットでちょっと調べてみたところ,鳩の糞はどうやら幸運のサインみたいです(笑)ということで,2021年は自分にとってかつてない豊作の年に,,,なればいいなあ.

やりたいことリスト!

ブログ書き始めて4年くらいになってまして,過去のエントリを振り返ってみてたんですが,「やりたいこと!」って書いてあることほとんど終わってないんですよね(笑)実際に結構時間かけてやってることもあるんですが,なかなかまとまったとこまでたどり着くのは大変で...なので,期限を切らずに(笑)じっくりやっていきたいと思います.

【仕事ネタ】

1.1人オフィスの近くに引っ越し.(流山市ー>市川市
 通勤時間がもったいないので,南流山から引っ越します!南流山自体は大分気に入ってたので残念なのですが.横浜に戻るのもありかと思ったんですが,都内への移動が時間かかるのでおそらく行徳とか浦安とかになりそうです.
2.会社のロゴ作成
 引っ越したら名前をちょっとづつ出していきます!
3.税理士を見つける.
4.初決算!
 初決算,自分でやってみようか悩んでたんですが,法人決算は結構難しいらしく,ミスって追徴課税におびえるのもつらいので,税理士を探すことにしました.ただ,キャッシュフローがめっちゃシンプルなので,ちょっともったいないなあ...という思いもあるのですが.
5.お客さん&売上を増やす!
 完全に仕事漬けの2021年にするか,,,どうするか...
6.海外から仕事をとってくる!
 今すぐは無理ですが,今後の目標に.

【勉強ネタ】

1.つくばチャレンジのスポンサーになる.
2.つくばチャレンジ完走.
 コロナ次第だと思うのですが,今年実施されれば参加します!今年は会社で参加する(=宣伝目的含む)ので,本気です(笑)
3.SFMをやりきる.
 うーん.やっている.が,終わらない.
4.ディープラーニングを勉強する.
 うーん.やりたい.が,始められない.
5.数学力強化!
 現在進行形です!近いうちに何かまとめを残したい(/・ω・)/
6.月に1本のブログアップ
 記憶の定着のために,,,,月に一つは勉強したことまとめたい.

でやりたいこと,見ての通り【仕事ネタ】と【勉強ネタ】に分かれるんですが,「1人会社でどこまで売り上げられるか!」に挑戦するか,「仕事をある程度の量に抑えて,勉強しつつ会社の名前を売っていく方向にもっていく」か...をちょっと悩んでいます.売上アゲアゲ方向にもっていくと,完全に1年間仕事漬けになることと,どうしても時間の切り売り的側面が強くなってしまう一方で,せっかく会社作ったので,どこまでお金稼げるかやってみたい気持ちもあり(笑)ただ,一回やり始めると途中下車できないので,悩み中です.

2021年!

ということで,今年もいままで通り,

All life is an experiment. The more experiments you make the better.
Ralph Waldo Emerson

をモットーに,いろんな実験をしていけたらと思います!ではでは,今年もよろしくお願いします!

f:id:rkoichi2001:20210102234243j:plain
弊社にさらなる幸運が舞い込みますように(/・ω・)/

ビックマック

前を向きたいとき,僕はビックマックを食べる.
仕事が終わらず,2.5㎡のオフィスでコードとにらめっこしていた時,僕はビックマックを食べた.
クリスマスに誰も出迎えてくれなかったときも,僕はビックマックを食べた.
「年末,休めないんじゃね..」と不安になりながらも,僕はビックマックを食べた.
この先も,何度もくじけそうになるだろう.だから,これからもよろしく.

f:id:rkoichi2001:20201230171635j:plain
僕はビックマックを食べる.


ということで,仕事終わったーーー!!!なんとか冬休みを迎えられそうです.
仕事納めのご褒美として,久々にビックマックを食べてみました.「ご褒美にマックって,どんな貧乏社長やねん」って言われそうですが,お金的にというよりは,カロリー&背徳感的なご褒美ですかね.

結局今年はあんまりブログ書けませんでしたが,,,,冬休みはちょっと頑張ってアップできればと思います(゚Д゚)ノ

オフィス

どうもどうも,社長です(笑).

近況

 ちょっと更新が滞ってしまったんですが,実は7月頭から会社が始動しておりまして,せっせと働き始めてます.だいぶスタートダッシュをかまさないといけない状況でありまして(笑),ブログの更新頻度がちょっと落ちてしまうと思うのですが,ちょこちょことアップしますのでまた暇なときに見ていただければありがたいですm(__)m
(今年はつくばチャレンジも不参加なので,SFMもなんとかやり切ります!)

オフィス契約

 勢いで法人を立ち上げ,一人会社なので自宅を登記したのですが,これ,国税庁のホームページで住所出てくるんですよね...登記後にそれに気づいたんですが,後の祭りでして...会社のホームページ立ち上げたり,ブログで会社名を公表しようと思ってたんですが,身に覚えのないピザとかがいっぱい届いたら困るので(笑)オフィスを借りて,住所を登記しなおすことにしました.

 で,まともなオフィスはさすがに高いので,レンタルオフィスを借りることにしました.残念ながら南流山近辺にはなかったので,ビズサークルという会社が運営している市川妙典のレンタルオフィスにしました.南流山からは少し離れているので,マンション更新のタイミングで市川市に引っ越そうと思います.ということで,なかなか港区民への道のりは遠そうです(笑)

 肝心のオフィスですが,共益費込みで月額 32000 円でして,何とかやりくりできそうな額です.下にオフィスの写真を載せてます.自分でいうのもなんですが,あまりの ”おままごと” っぷりに自然と笑みがこぼれます(爆)が,今をときめく Apple もガレージから始まったことを考えると,先はどうなるかわかりませんよ~.目指せ100兆円起業 ( ゚Д゚)

f:id:rkoichi2001:20200718223403j:plain
市川妙典の Anytime Fitness! ではなくて,オフィスです.
f:id:rkoichi2001:20200718223500j:plain
映画マトリックスを彷彿とさせる,ドアの群れ.
f:id:rkoichi2001:20200718223602j:plain
わが社のオフィス!ここで缶詰になってます

訪問調査

 ある程度の規模の会社になると,初回取引では相手法人の信用調査をします.ただ,法人調査部門を会社毎に持つことは難しいので,普通は調査会社に外注します.日本の場合だと,東京商工リサーチ帝国データバンクの2社が有名どころだと思います.
 実は,会社(=自宅)に東京商工リサーチの調査員の方が来ました!訪問調査というよりは,家庭訪問みたいな形になってしまったんですが,一応となりのコンビニでおーいお茶を購入して,丁重にもてなしました.ということで僕の会社,企業情報データベースに乗ってます!スコアはおそらく最低レベルだと思いますが(笑)

それでは,今後ともわが社をよろしくお願いしますm(__)m
(まだ先になりそうですが,会社住所の登記変更が完了したら,社名&ホームページも公開します!)

Fisher Vector を用いた類似画像検索

ということで,ちょっと時間がかかってますが類似画像検索まで来ました!

X. 類似画像検索

類似画像検索とは何か?という話なんですが,Googleの画像検索のように,クエリ画像に対して「よく似ている」画像を検索することです.

「類似画像検索がなぜSFMに必要なのか?」という話なんですが,SFMを実行する場合,ペアとなる写真には同じ物体が写っている必要があります.例えば東京タワーを3次元復元したい場合,東京タワーが写っていない写真を使っても意味がありません.もし手持ちの写真に東京タワーが写っていない写真も含まれているとしたら,これを除く必要があります.下記の例では,スカイツリーの写真を除外しています.

f:id:rkoichi2001:20200617222657p:plain
東京タワーとスカイツリー

復元したい物体のスケールが「東京タワー」とか「奈良の大仏」程度なら,頑張れば一枚の写真に全てを収めることができます.が,「ローマの町並み」とか「国際通り」とかになってくると,写真一枚に全てを収めることは当然できないので,オーバーラップのある写真をペアとして繋いでいく必要があります.

f:id:rkoichi2001:20200617224702p:plain
国際通りのお菓子御殿

上の写真は国際通りのお菓子御殿(笑)ですが,3つの写真を使ってお菓子御殿を復元したい場合,左の写真は少なくとも中央の写真と,右の写真も少なくとも中央の写真とペアにならなければ3枚の写真間につながりができません.ということで,類似している画像を検索することによって,クエリ画像とペアになる写真を見つけていきます.(※写真が時系列に並んでいる場合,おそらく直近の数枚の写真がペアになるので,類似検索は必ずしも必要無いです.)

X. Fisher Vector とは

いつもの如く TheiaSfM で使われていた方法を参考に, Fisher Vector による方法を用いました.Fisher Vector の原理を理解するには,統計とか情報幾何学とかをかじる必要があって,数学力の低い私は今頭を抱えているのですが...このへんもある程度理解できたらまたエントリを残せればと思います.

原理は難しいものの,ライブラリを使えば処理フローはそれほど複雑ではなかったので,簡単に書き留めておきます.大まかに4つの処理からなっていまして,

1.全画像の全SIFT特徴量を用いて,GMM(混合ガウスモデル)を生成.
 ー>生成過程では,全画像のSIFT特徴量を入力してやればOKです.多すぎる場合,間引いた方がいいかもしれません.この混合ガウスモデルが特徴量の生成モデルになります.

f:id:rkoichi2001:20200618000652j:plain
生成モデルを作成

2.生成した GMM を用いて,個々の画像を Fisher Vector に変換.
 ー>一つ一つの画像データを Fisher Vector に変換します.ここではFisher Vector が一つの画像を表現しています.今回の設定では一つの画像から大体15000個位の SIFT 特徴量が生成されたのですが,各識別子の次元が 128 なので,SIFT を用いて一つの画像を表現すると,15000 x 128 = 1920000 位の次元のベクトルになります.が,これでは大きいので,生成した GMM のパラメータ(の勾配)空間を使って次元削減してやります.この次元削減によって,今回のケースだと一枚の画像が 8000 次元くらいのベクトル(= Fisher Vector)で表現できるようになりました.

f:id:rkoichi2001:20200618000737j:plain
画像データを Fisher Vector に変換

3.全画像ペア間の Fisher Vector の距離を計算します.
 ー>「2つの画像が似ている」=「2つの画像の Fisher Vector が似ている」ということになります.あとは,全ペアの Fisher Vector の距離を総当り計算して行列として保持しておけば,画像間の類似関係をいつでも参照することができます.

f:id:rkoichi2001:20200618000803j:plain
Fisher Vector 間(=画像間)の距離を Matching Matrix として保持.

4.クエリ画像とその他画像の距離チェックし,距離が近いものを類似画像とする.
 ー> クエリ画像のインデックスをもつ Matching Matrix の列(もしくは行)を取り出せば,クエリ画像とその他全画像の距離がわかります.

X. 生成した SIFT 識別子による GMM のトレーニン

下記,ほぼ VLFeat ライブラリの呼び出しをしただけですが,,,

void FisherVectorEncoder::TrainGMM(uint64_t max_itr,
                                   const std::vector<Eigen::VectorXf>& train_data) {
  // X. Create GMM object.
  LOG(INFO) << "vl_gmm_new : ";
  VlGMM* gmm = vl_gmm_new(VL_TYPE_FLOAT, num_dimension_, num_clusters_);

  // X. Set the maximum number of EM iterations to 100.
  LOG(INFO) << "vl_gmm_set_max_num_iterations : ";
  vl_gmm_set_max_num_iterations(gmm, max_itr);

  // X. Set the initialization to random selection.
  LOG(INFO) << "vl_gmm_set_initialization : ";
  vl_gmm_set_initialization(gmm, VlGMMRand);

  // X. Cluster the data. (Train the GMM).
  num_data_ = train_data.size();
  LOG(INFO) << "Converting to MatrixXf : ";
  Eigen::MatrixXf mat_data = ConvertVectorOfFeaturesToMatrix(train_data);

  LOG(INFO) << "Data count for training is : " << num_data_;
  vl_gmm_cluster(gmm, mat_data.data(), num_data_);

  LOG(INFO) << "Training is Done! Now storing data!";
  LOG(INFO) << "means copied!";
  means_.resize(num_dimension_ * num_clusters_);
  std::memcpy(means_.data(), vl_gmm_get_means(gmm), sizeof(float) * num_dimension_ * num_clusters_);

  LOG(INFO) << "covariances copied!";
  covariances_.resize(num_dimension_ * num_clusters_);
  std::memcpy(covariances_.data(), vl_gmm_get_covariances(gmm),
              sizeof(float) * num_dimension_ * num_clusters_);

  LOG(INFO) << "priors copied!";
  priors_.resize(num_clusters_);
  std::memcpy(priors_.data(), vl_gmm_get_priors(gmm), sizeof(float) * num_clusters_);

  LOG(INFO) << "posteriors copied!";
  posteriors_.resize(num_clusters_ * num_data_);
  std::memcpy(posteriors_.data(), vl_gmm_get_posteriors(gmm),
              sizeof(float) * num_clusters_ * num_data_);

  LOG(INFO) << "Deleting gmm";
  vl_gmm_delete(gmm);
  LOG(INFO) << "gmm deleted.";
}

X. 画像間距離マトリクスの作成

全画像間の距離を計算しています.ここも結構時間がかかります.

void CreateMatchingMatrix(const std::string& gmm_file_path, const std::string& matching_matrix_path,
                          std::unordered_map<std::string, ImageInfo>& image_info_map) {
  FisherVectorEncoder fisher_encoder(gmm_file_path);

  // X.
  LOG(INFO) << "Create filename index map.";
  std::unordered_map<std::string, int> filenames_indices;
  std::unordered_map<int, std::string> indices_filemanes;
  int index = 0;
  for (auto citr = image_info_map.cbegin(); citr != image_info_map.cend(); ++citr) {
    std::filesystem::path p(citr->first);
    std::string filename = p.replace_extension("").filename().string();
    filenames_indices[filename] = index;
    indices_filemanes[index] = filename;
    index++;
  }

  LOG(INFO) << "Matching matrix computation starts!";
  int size = image_info_map.size();
  Eigen::MatrixXf matching_matrix(size, size);

  int cnt = 0;
  int total_size = image_info_map.size();
  for (const auto& [filename1, value] : image_info_map) {
    int idx1 = filenames_indices[filename1];
    matching_matrix(idx1, idx1) = 0.0;

    LOG(INFO) << "Image Distance... " << cnt << " / " << total_size;
    cnt++;
    for (int idx2 = idx1 + 1; idx2 < size; idx2++) {
      Eigen::VectorXf feature1 =
          fisher_encoder.ComputeFisherVector(image_info_map[filename1].descriptors_);

      std::string filename2 = indices_filemanes[idx2];
      Eigen::VectorXf feature2 =
          fisher_encoder.ComputeFisherVector(image_info_map[filename2].descriptors_);

      float score = (feature1 - feature2).squaredNorm();

      matching_matrix(idx1, idx2) = score;
      matching_matrix(idx2, idx1) = score;
    }
  }

  std::cout << matching_matrix << std::endl;

  std::ofstream writer(matching_matrix_path, std::ios::out | std::ios::binary);
  cereal::PortableBinaryOutputArchive output_archive(writer);
  output_archive(filenames_indices, matching_matrix);
}

X. 実行結果!

ということで,類似画像検索をやってみました!使った画像セットは,SIFTのエントリのときにリンクを貼った下記の画像セットです.

https://research.cs.cornell.edu/1dsfm/
上記サイトの Gendarmenmarkt images (tar, 1.0 GB) をダウンロードして使いました.

クエリ画像としては,下記の画像を使いました.

f:id:rkoichi2001:20200619110826j:plain
クエリ画像

で結果ですが,スコアの低い(類似度が高い)画像から順に上位10枚の画像が下記のようになりました.下記の結果を見てもわかるのですが,やっぱり天候とか昼夜の条件によってスコアに開きが出てしまうのは仕方ないのかもしれませんね...この辺の影響を抑える方法があれば,つくばチャレンジにも使えるようになるかも.

f:id:rkoichi2001:20200619113439p:plain
クエリ画像と上位10画像

スコアの分布は下記のようになります.最低スコア(類似度が高い)0.516687を皮切りに,上限2.3くらいまで連続的に分布していきます.

f:id:rkoichi2001:20200619111602p:plain
スコア分布

X. 実験に使ったコード

github.com

VLFeat & Cereal も Submodule として取り込んであるので,下記のコマンドを叩いて貰えれば動くとこまでは行きます.

git clone --recurse-submodules git@github.com:koichiHik/blog_vlfeat_image_compare.git
cd blog_vlfeat_image_compare/
sh build.sh

で,今回は実行ファイルが3つありまして,順を追って実行する必要があります.
1.作成済のSIFT特徴量を用いて,GMMの生成
(flagfile.txt に SIFT 特徴量,画像ファイルのディレクトリパスを指定し,生成した GMM モデルのパラメータファイルの出力パスを指定します.)

./build/train_gmm_model --flagfile=./flagfiles/train_gmm_model.txt

2.生成したGMMを用いて,全画像の Fisher Vector を計算.各画像間の距離マトリクスを作成する.
(flagfile.txt に SIFT 特徴量のディレクトリパス,生成した GMM のパラメータファイルパスを指定し.生成するマッチングマトリクスの出力ファイルパスを指定します.)

./build/create_matching_matrix --flagfile=./flagfiles/train_gmm_model.txt

3.クエリ画像のファイル名を指定すると,類似している上位画像を出力します.
(flagfile.txt に SIFT 特徴量,画像ファイルのディレクトリパス,マッチングマトリクスのファイルパスと,ターゲット画像のファイル名を指定します.)

./build/show_similar_images --flagfile=./flagfiles/show_similar_images.txt

X. 参考文献

- Exploiting generative models in discriminative classifiers. 

- Fisher Kernels on Visual Vocabularies for Image Categorization. 

- Image Classification with the Fisher Vector : Theory and Practice. 

VLFeat ライブラリを使った SIFT の計算

またまた SIFT のエントリになってしまいますが,TheiaSfM の特徴点計算で OpenCV ではなくて VLFeat が使われていたこともあって,試しに使ってみることにしました.どこかの論文に書いてあったと思うのですが,OpenCV の SIFT よりも VLFeat の SIFT のほうが特徴点がたくさんでるとか出ないとか...ここで生成した SIFT を保存して,いよいよ次のエントリで画像間の類似度計算をします.

X. VLFeat

VLFeat ライブラリは画像識別や特徴点の抽出/マッチングに特化したアルゴリズムのライブラリになります.VLFeat で実装されているアルゴリズムの中には,OpenCV に実装の無いものもあって,この後の画像類似度計算のところで使う Fisher Vector 生成の実装も現状では OpenCV には無いような気がします.ドキュメント周りも結構整備してくれてありまして,例えば SIFT の計算のチュートリアルは下記になります.

www.vlfeat.org

X. VLFeat の SIFT の使い方

で,実際に VLFeat を使った SIFT 特徴量の計算ですが,大きな流れは下記のようになります.次の Octave に行くかどうかがライブラリユーザに委ねられてるところとか,ちょっと面白い作りだなと思いました.

0. u8 型から float 型に型変更.
  // 0. Convert image to float.
  cv::Mat float_img(gray_image.rows, gray_image.cols, CV_32FC1);
  gray_image.convertTo(float_img, CV_32FC1);
1. 1st オクターブをどうするか決定します.この値が負の値なら,入力画像をアップサンプル(拡大)します.
  // 1. Compute Valid First Octave.
  const int first_octave = GetValidFirstOctave(params.def_first_octave_, float_img.cols,
                                               float_img.rows, params.max_scaled_dim_);
2. VLFeat の SIFT 生成器を生成します.(vl_sift_new)
  // 2. Create new sift filter.
  VlSiftFilt* sift_filter = vl_sift_new(float_img.cols, float_img.rows, params.num_octave_,
                                        params.num_levels_, first_octave);
3. 必要なパラメータを設定.
  // 3. Set Parameters.
  vl_sift_set_peak_thresh(sift_filter, params.peak_threshold_);
  vl_sift_set_edge_thresh(sift_filter, params.edge_threshold_);
  vl_sift_set_norm_thresh(sift_filter, params.norm_threshold_);
  vl_sift_set_magnif(sift_filter, params.magnif_);
  vl_sift_set_window_size(sift_filter, params.window_size_);
4. 1st オクターブを処理します.
  // 4. Process First Octave.
  int octave_no = 1;
  int vl_status =
      vl_sift_process_first_octave(sift_filter, reinterpret_cast<float*>(float_img.data));
5.1. 現在のオクターブに関して,キーポイントを計算・取得します.
    // Detect Key Points.
    vl_sift_detect(sift_filter);
5.2. 取得したキーポイント一つ一つに対してループをまわし,識別子を計算する.
      // Compute orientations.
      double angles[4];
      int num_angles = vl_sift_calc_keypoint_orientations(sift_filter, angles, &vl_keypoints[i]);

      // Compute descriptor for all keypoints and orientations.
      Eigen::VectorXf descriptor(params.num_sift_dimensions_);
      for (int j = 0; j < num_angles; ++j) {
        descriptor.setZero();
        vl_sift_calc_keypoint_descriptor(sift_filter, descriptor.data(), &vl_keypoints[i],
                                         angles[j]);
        descriptors.push_back(descriptor);
        KeyPoint kp(vl_keypoints[i].x, vl_keypoints[i].y, vl_keypoints[i].sigma, angles[j], 0,
                    vl_keypoints[i].o);
        keypoints.push_back(kp);
      }
5.3. 次のオクターブへ向かう.
    // Proceed to next octave.
    vl_status = vl_sift_process_next_octave(sift_filter);
    octave_no = octave_no + 1;
6. SIFT 生成器を削除し,終了.
  // X. Delete sift caculator.
  vl_sift_delete(sift_filter);

X. VLFeat の SIFT の使い方(実際の実装)


void ComputeKeyPointsAndDescriptors(const cv::Mat& gray_image, const SIFTParams& params,
                                    std::vector<KeyPoint>& keypoints,
                                    std::vector<Eigen::VectorXf>& descriptors) {
  keypoints.clear();
  descriptors.clear();

  // 0. Convert image to float.
  cv::Mat float_img(gray_image.rows, gray_image.cols, CV_32FC1);
  gray_image.convertTo(float_img, CV_32FC1);

  // 1. Compute Valid First Octave.
  const int first_octave = GetValidFirstOctave(params.def_first_octave_, float_img.cols,
                                               float_img.rows, params.max_scaled_dim_);

  // 2. Create new sift filter.
  VlSiftFilt* sift_filter = vl_sift_new(float_img.cols, float_img.rows, params.num_octave_,
                                        params.num_levels_, first_octave);

  // 3. Set Parameters.
  vl_sift_set_peak_thresh(sift_filter, params.peak_threshold_);
  vl_sift_set_edge_thresh(sift_filter, params.edge_threshold_);
  vl_sift_set_norm_thresh(sift_filter, params.norm_threshold_);
  vl_sift_set_magnif(sift_filter, params.magnif_);
  vl_sift_set_window_size(sift_filter, params.window_size_);

  // 4. Process First Octave.
  int octave_no = 1;
  int vl_status =
      vl_sift_process_first_octave(sift_filter, reinterpret_cast<float*>(float_img.data));

  // 5. Process Octave until we cant anymore.
  while (vl_status != VL_ERR_EOF) {
    // Detect Key Points.
    vl_sift_detect(sift_filter);

    // Get Key Points.
    const VlSiftKeypoint* vl_keypoints = vl_sift_get_keypoints(sift_filter);
    const int num_keypoints = vl_sift_get_nkeypoints(sift_filter);

    // Compute descriptors for all keypoints detected.
    for (int i = 0; i < num_keypoints; ++i) {

      // Compute orientations.
      double angles[4];
      int num_angles = vl_sift_calc_keypoint_orientations(sift_filter, angles, &vl_keypoints[i]);

      // Compute descriptor for all keypoints and orientations.
      Eigen::VectorXf descriptor(params.num_sift_dimensions_);
      for (int j = 0; j < num_angles; ++j) {
        descriptor.setZero();
        vl_sift_calc_keypoint_descriptor(sift_filter, descriptor.data(), &vl_keypoints[i],
                                         angles[j]);
        descriptors.push_back(descriptor);
        KeyPoint kp(vl_keypoints[i].x, vl_keypoints[i].y, vl_keypoints[i].sigma, angles[j], 0,
                    vl_keypoints[i].o);
        keypoints.push_back(kp);
      }
    }

    // Proceed to next octave.
    vl_status = vl_sift_process_next_octave(sift_filter);
    octave_no = octave_no + 1;
  }

  // X. Delete sift caculator.
  vl_sift_delete(sift_filter);
}

X. この後....

このコードで一応画像毎に SIFT のキーポイントと識別子のファイル出力ができるようになってるんですが,ここで生成した SIFT を使って次のエントリで画像類似度を計算します.使った画像はCornell 大の下記のページからダウンロードできます.
https://research.cs.cornell.edu/1dsfm/

上記サイトの Gendarmenmarkt images (tar, 1.0 GB) をダウンロードして使いました.

X. 実験に使ったコード

github.com

VLFeat & Cereal も Submodule として取り込んであるので,下記のコマンドを叩いて貰えれば動くとこまでは行きます.

git clone --recurse-submodules  git@github.com:koichiHik/blog_vlfeat_sift.git
sh build.sh
./build/sift_extractor --flagfile=./flagfiles/sift_test_flags.txt

X. 参考文献&プロジェクト

TheiaSfM
github.com

vlfeat
github.com

特徴点(SIFT)のシリアライズとその他もろもろ

画像間類似度を計算する方法を調査してるんですが,画像の枚数が多くなってくると毎回SIFTを計算するのも辛くなってくるので,結果を保存する方法を調べることにしました.このトピックに限らずですが,参考にしたのはTheiaSfMというオープンソースプロジェクトです.githubのリンクは最後に載せてます.

X.Cereal

boost とか,いくつかシリアライズのためのライブラリがあるみたいなんですが TheiaSfM で Cereal が使われていたこともあって,ひとまず使ってみることにしました.ドキュメントも用意してくれていて,とても親切な感じですm(_ _)m

簡単なチュートリアル
cereal Docs - Main

githubのリンク
github.com

X.今回のシリアライズ対象.

で,今回シリアライズしてファイル保存したい内容ですが,SIFTのキーポイントと識別子になります.キーポイントの定義は,,,

class KeyPoint {
 public:
  KeyPoint() {}
  KeyPoint(float x, float y, float size, float angle, float responce, int octave)
      : x_(x), y_(y), size_(size), angle_(angle), responce_(responce), octave_(octave) {}

  bool operator==(const KeyPoint& rhs) const {
    return x_ == rhs.x_ && y_ == rhs.y_ && size_ == rhs.size_ && angle_ == rhs.angle_ &&
           responce_ == rhs.responce_ && octave_ == rhs.octave_;
  }

  bool operator!=(const KeyPoint& rhs) const { return !this->operator==(rhs); }

  friend std::ostream& operator<<(std::ostream& output, const KeyPoint& kp) {
    output << kp.x_ << ", " << kp.y_ << ", " << kp.size_ << ", " << kp.angle_ << ", "
           << kp.responce_ << ", " << kp.octave_;
    return output;
  }

 private:
  float x_;
  float y_;
  float size_;
  float angle_;
  float responce_;
  int octave_;

  friend cereal::access;
  template <class Archive>
  void serialize(Archive& ar) {
    ar(x_, y_, size_, angle_, responce_, octave_);
  }
};

で,識別子の定義は Eigen と STLVector を使って表現します.

std::vector<Eigen::VectorXf>& descriptors

X.Cerial をつかって自前のクラス(KeyPoint)を保存する.

上記で定義した KeyPoint クラスを Cereal を使って Serialize する方法ですが,これはとても簡単で,上記クラス定義の中の下記の関数を記述してあげればOKです.

  friend cereal::access;
  template <class Archive>
  void serialize(Archive& ar) {
    ar(x_, y_, size_, angle_, responce_, octave_);
  }

X.Cerial をつかって Eigen を保存する.

こちらがちょっとむずかしくて理解するのに時間がかかったんですが,シリアライズしたいオブジェクトを引数に取る "load" と "save" という関数を定義します.

template <class Archive, class _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows,
          int _MaxCols>
inline typename std::enable_if<traits::is_output_serializable<BinaryData<_Scalar>, Archive>::value,
                               void>::type
save(Archive& ar, Eigen::Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols> const& m) {
  int32_t rows = m.rows();
  int32_t cols = m.cols();
  ar(rows);
  ar(cols);
  ar(binary_data(m.data(), rows * cols * sizeof(_Scalar)));
}

template <class Archive, class _Scalar, int _Rows, int _Cols, int _Options, int _MaxRows,
          int _MaxCols>
inline typename std::enable_if<traits::is_input_serializable<BinaryData<_Scalar>, Archive>::value,
                               void>::type
load(Archive& ar, Eigen::Matrix<_Scalar, _Rows, _Cols, _Options, _MaxRows, _MaxCols>& m) {
  int32_t rows;
  int32_t cols;
  ar(rows);
  ar(cols);
  m.resize(rows, cols);
  ar(binary_data(m.data(), static_cast<std::size_t>(rows * cols * sizeof(_Scalar))));
}

上記のコードは TheiaSfM から持ってきたんですが,Eigen のコンストラクタが取るテンプレートパラメータを "load" と "save" に全て引き継がせることで,任意の Eigen::Matrix をシリアライズ/デシリアライズできるようにしてます.で,なおかつ std::enable_if を用いることで,Eigen::Matrix の中身の要素が Ceral によってバイナリ化できるかどうか調べてます.なので,,,バイナリ化できなければコンパイル時に関数が定義されず,エラーになる...はず...

C++得意な人って,このあたりって「息をする」用にコーディングできるんですかね....牛歩プログラマとしては羨ましい限りです.

X.実験に使ったコード

Cereal の github も submodule として取り込んでいるので,下記のコマンドを叩いてもらえば動くとこまでは行きます.

git clone --recurse-submodules  git@github.com:koichiHik/blog_cereal_use.git
sh build.sh
./build/cereal_sample 

github.com

X.参考文献&プロジェクト

TheiSfM
github.com

X.実際のコード

で,下記が実験用のコードです.単純にシリアライズ/デシリアライズして入出力前後のデータが一致していることを見てるだけですが,,,,


// STL
#include <cstdio>
#include <fstream>
#include <iostream>
#include <ostream>
#include <string>
#include <vector>

// Cereal
#include <cereal/archives/portable_binary.hpp>
#include <cereal/types/vector.hpp>

// Glog
#include <glog/logging.h>

// Eigen
#include <eigen3/Eigen/Core>

// OpenCV
#include <opencv2/features2d/features2d.hpp>

// Original
#include "eigen_serializable.h"

class KeyPoint {
 public:
  KeyPoint() {}
  KeyPoint(float x, float y, float size, float angle, float responce, int octave)
      : x_(x), y_(y), size_(size), angle_(angle), responce_(responce), octave_(octave) {}

  bool operator==(const KeyPoint& rhs) const {
    return x_ == rhs.x_ && y_ == rhs.y_ && size_ == rhs.size_ && angle_ == rhs.angle_ &&
           responce_ == rhs.responce_ && octave_ == rhs.octave_;
  }

  bool operator!=(const KeyPoint& rhs) const { return !this->operator==(rhs); }

  friend std::ostream& operator<<(std::ostream& output, const KeyPoint& kp) {
    output << kp.x_ << ", " << kp.y_ << ", " << kp.size_ << ", " << kp.angle_ << ", "
           << kp.responce_ << ", " << kp.octave_;
    return output;
  }

 private:
  float x_;
  float y_;
  float size_;
  float angle_;
  float responce_;
  int octave_;

  friend cereal::access;
  template <class Archive>
  void serialize(Archive& ar) {
    ar(x_, y_, size_, angle_, responce_, octave_);
  }
};

void CreateKeyPointAndDescriptorPair(int num_data, std::vector<KeyPoint>& keypoints,
                                     std::vector<Eigen::VectorXf>& descriptors) {
  keypoints.clear();
  descriptors.clear();

  for (int i = 0; i < num_data; i++) {
    KeyPoint kp(1.1 * i, 1.1 * 10 * i, 1.1 * 100 * i, 1.1 * 1000 * i, 1.1 * 10000 * i,
                1.1 * 100000 * i);
    Eigen::VectorXf vec(4);
    vec << 1.1 * i, 1.1 * 10 * i, 1.1 * 100 * i, 1.1 * 1000 * i;
    keypoints.push_back(kp);
    descriptors.push_back(vec);
  }
}

void SaveAsBinary(const std::string& filepath, const std::vector<KeyPoint>& keypoints,
                  const std::vector<Eigen::VectorXf>& descriptors) {
  std::ofstream feature_writer(filepath, std::ios::out | std::ios::binary);
  cereal::PortableBinaryOutputArchive output_archive(feature_writer);
  output_archive(keypoints, descriptors);
}

void LoadFromBinary(const std::string& filepath, std::vector<KeyPoint>& keypoints,
                    std::vector<Eigen::VectorXf>& descriptors) {
  keypoints.clear();
  descriptors.clear();

  std::ifstream feature_reader(filepath, std::ios::in | std::ios::binary);
  cereal::PortableBinaryInputArchive input_archive(feature_reader);

  input_archive(keypoints, descriptors);
}

bool Compare(const std::vector<KeyPoint>& org_keypoints,
             const std::vector<Eigen::VectorXf>& org_descriptors,
             const std::vector<KeyPoint>& loaded_keypoints,
             const std::vector<Eigen::VectorXf>& loaded_descriptors) {
  for (int i = 0; i < org_keypoints.size(); i++) {
    if (org_keypoints[i] != loaded_keypoints[i] || org_descriptors[i] != loaded_descriptors[i]) {
      return false;
    }
    std::cout << std::endl;
    std::cout << org_keypoints[i] << std::endl;
    std::cout << loaded_keypoints[i] << std::endl;
    std::cout << org_descriptors[i] << std::endl;
    std::cout << loaded_descriptors[i] << std::endl;
  }
  return true;
}

int main(int argc, char** argv) {
  google::ParseCommandLineFlags(&argc, &argv, true);
  FLAGS_alsologtostderr = 1;
  FLAGS_stderrthreshold = google::GLOG_INFO;
  google::InitGoogleLogging(argv[0]);

  // X. Setting.
  const std::string test_path = "./sample.bin";
  int num_data = 10;
  std::vector<KeyPoint> keypoints;
  std::vector<Eigen::VectorXf> descriptors;

  // X. Create keypoints and descriptors.
  CreateKeyPointAndDescriptorPair(10, keypoints, descriptors);

  // X. Serialize and save to file.
  SaveAsBinary(test_path, keypoints, descriptors);

  // X. Load from file and deserialize.
  std::vector<KeyPoint> loaded_keypoints;
  std::vector<Eigen::VectorXf> loaded_descriptors;
  LoadFromBinary(test_path, loaded_keypoints, loaded_descriptors);

  // X. Compare if the results are same.
  if (!Compare(keypoints, descriptors, loaded_keypoints, loaded_descriptors)) {
    LOG(INFO) << "Contents are not the same!";
  } else {
    LOG(INFO) << "Contents are the same!";
  }

  // X. Remove temporary file.
  std::remove(test_path.c_str());

  return 0;
}

最近傍探索(Nearest Neighbor Search)

ということで,なんだかんだ会社設立やその他もろもろが思った以上に大変で,,,ちょっと止まってたんですが,ひとまず小休止に入ったので引き続き SFM をやりたいと思います!

X. 最近傍探索とは.

まず,最近傍探索とは何ぞや?って話なんですが,Wikipedia の定義を見ると,
距離空間における最も近い点を探す最適化問題の一種」
とあります.「距離空間」には「マンハッタン距離」とか「ユークリッド距離」とかいろいろとありますが,ここでは簡単のために「ユークリッド距離」前提としてエントリを進めていきます.で,これだけ見ても意味わからないので,まず簡単に一次元の例を考えてみたいと思います.

Brute Force Macher 1次元データの最近傍探索

f:id:rkoichi2001:20200606200013j:plain
力任せ法による最近傍探索

この問題は,「クエリ点(=4.2)にもっとも近い値を探索対象の中から探す.」ということになりますが,上記は力任せ方による解法です.このやり方だとクエリ点と探索対象点のすべての距離を計算して,距離が最も近いものを選択するという流れになります.このやり方では,比較対象点がn個存在する場合,探索回数が O(n) になります.なので,比較対象点が1億点あるときは1億回の比較が必要になります.

Binary Treeによる1次元データの最近傍探索

f:id:rkoichi2001:20200606200051j:plain
Binary Tree による最近傍探索

次に,Binary Tree を用いた最近傍探索になります.力任せ探索とは異なり,比較対象点が n 個存在する場合,探索回数が O(log n) のオーダになります.実際には数値のばらつき具合によりますが,うまく行けば比較対象点が1億点あるときは,約27回程度の比較になります.めちゃ減りましたね.上図では,ツリーの葉ノードが探索対象を表しており,四角のノードが分岐のための閾値を表しています.閾値は単純に当該ノード以下の探索対象の平均値を用いています.

で,今回のエントリでは,「点が1次元じゃなくて,高次元のユークリッド空間にある場合はどうやって探索すんの?」という話になります.高次元になったときでも力任せ探索は有効で,クエリ点に対して一番近い点を確実に見つけてくれます.ただ,時間がかかります.一方で, Binary Tree を高次元に適用しようとすると,ちょっと工夫が必要そうです.

X. SFMの処理フローと計算コスト.

で,過去に書いた SFM エントリのちょっとしたおさらいになるんですが,細かな処理を除くと SFM の大きな処理フローとは下記になります.

1.各写真から特徴点の抽出.

f:id:rkoichi2001:20200606135642p:plain
特徴点の抽出

2.視点が似ている(=同じようなものが被写体になっている)写真のペアを見つける.
【探索コスト大】全ての写真ペアを比較し,視点が似ているかチェック.

f:id:rkoichi2001:20200606140853p:plain
写真ペアの探索

3.見つけた写真のペアに紐づく特徴点をマッチング.
【探索コスト大】ペアの全ての特徴点を比較して,最近傍ペアを推定する.

f:id:rkoichi2001:20190721222640p:plain
特徴点ペアの探索

4.カメラの位置と点群を推定する.

f:id:rkoichi2001:20190721230139p:plain
カメラ位置と点群の推定

で,2と3のところで探索の話が出てきます.2では写真ペア,3では特徴点ペアということで,ちょっと探索の内容は異なりますが,SFMの復元対象物体が大きくなるほどこの2と3にかかるコストが大きくなっていきます.奈良の大仏とかその位のスケールの構造物でも写真の枚数がある程度大きくなるはずなので,2と3を工夫しないと計算量が爆発します.今回のエントリでは,3の特徴点探索の話をします.

※ちなみに那覇国際通りは1.6キロあるので,2と3を工夫するだけだとまだ不十分で,Hybrid SfM か,写真の時系列順序を使う方法を試したいと思ってます.

X. SIFTを用いた場合の SFM の特徴点マッチング.

今回は特徴点に SIFT を使おうと思っていますが, SIFT では画像中で見つかったキーポイントに一つづつに識別子(Descriptor)が付随しています.

f:id:rkoichi2001:20200606203252j:plain
SIFTのキーポイントと識別子

つまり,生成したキーポイントの比較に識別子を使うことになるのですが,この識別子は上図にあるように128個のビンを持つヒストグラムになっています.つまり,「128次元のユークリッド空間内に存在する点の最近傍探索」になります.

X. 今回解剖した3つの手法.

一つ一つのエントリがそこそこの分量になりそうだったので,別のポストにします!書き終わったら順次リンクを貼っていきます.

Brute Force Matching

エントリを書き終わったらリンクしますm(_ _;)m

KD Tree Search

エントリを書き終わったらリンクしますm(_ _;)m

Cascade Hashing for 3D Search

エントリを書き終わったらリンクしますm(_ _;)m

X. 3つの手法による最近傍探索の計算時間.

計算時間を全く考慮に入れなければ Brute Force Matching で良いのですが,そういうわけにもいかないので,計算時間とマッチング精度のトレードオフを取ることになります.で,今回は実験を簡単にするために,同じ画像から抽出したキーポイントと識別子を比較しました.ただ,比較対象となる特徴点セットの中に必ず完全一致する特徴点が見つかります.これでは実際に SFM において使う用途と状況が異なってしまうので,探索対象データに少しずつノイズを混ぜていって,探索時間や精度がどう変化するか見てみました.

実験に用いた画像

f:id:rkoichi2001:20200607192800p:plain
実験に用いた画像

取得できた特徴点数:2674

σ = 0 の場合の探索精度と時間
Brute Force Matcher ー> 処理時間:1.45s, 正解率:1.0
KDTree Matcher   ー> 処理時間:0.0062s, 正解率:1.0
Cascading Hasher  ー> 処理時間:0.129s, 正解率:1.0

σ = 4 の場合の探索精度と時間
Brute Force Matcher ー> 処理時間:1.45s, 正解率:1.0
KDTree Matcher   ー> 処理時間:0.0062s, 正解率:0.9996
Cascading Hasher  ー> 処理時間:0.129s, 正解率:0.9996

σ = 32 の場合の探索精度と時間
Brute Force Matcher ー> 処理時間:1.44s, 正解率:0.99
KDTree Matcher   ー> 処理時間:0.0252s, 正解率:0.75
Cascading Hasher  ー> 処理時間:0.1104s, 正解率:0.48

ということで,「ノイズが大きくなれば,KDTree Matcher の探索時間が大きくなり,Cascading Hasherのほうが性能がよくなります!!!」という結論に持って行きたかったんですが,どうやらそうならず(笑)ここはちょっと自分の理解不足な部分があるかもしれませんが,KDTree Matcher の性能が一番よい結果になりました.ただし,あくまでも同じ画像から生成した特徴点を比較しているので,この前提を変更するとひょっとすると結果が変わるかもしれません.

X. 実験に使ったコード

github.com

git clone して,下記のコマンドを実行すると上のテストが動きます.

sh build.sh
./build/nearest_neighbor_perf_test

X. 参考文献

KD-Trees に関連する論文.

An algorithm for Finding Best Matches in Logarithmic Expected Time.

Shape Indexing Using Approximate Nearest-Neighbor Search in High-Dimensional Spaces.

FAST APPROXIMATE NEAREST NEIGHBORS WITH AUTOMATIC ALGORITHM CONFIGURATION.

Optimized KD-trees for fast image descriptor matching.

Cascade Hashing に関連する論文.

Fast and Accurate Image Matching with Cascade Hashing for 3D Reconstruction.