0.バンドル調整 〜 目次
ご無沙汰しております!社長です!「月に1本のブログアップ」といった新年の誓いもなんのその(笑),更新せずに半年たっちゃいましたが,せっかくの自粛GWを活かして勉強したことまとめ書きします!ということで,以下,エントリです∠(`・ω・´)
春だ,桜だ,バンドル調整だ~.ということで,牛歩のように進めている一人SFMプロジェクト,やっとこさバンドル調整まで来ました(笑)仕事でも最適化が使いこなせるようになりたいと思っていたことも有り,少しずつ勉強していたのですが,ちょうどいい機会なのでバンドル調整の文脈で最適化のエントリを書いてみたいと思います.このページは目次として書いていますが,個々のページが書き終わった段階でここにリンクを貼りたいと思います.心が折れなければ,,,,最後まで書ききれると思います!
ちなみに,どの文献をよんでも「ヘッシアンの計算は大変なので,,,」とあったので,力ずくで計算してみました.どろんこになりながら実装したので,少しは最適化を身近に感じるようになりました(笑)
X.モデル化,及び問題設定
バンドル調整の簡単な説明,モデル化,最小化する数式を説明します.
daily-tech.hatenablog.com
X.回転のパラメータ化
回転行列のエントリは9つ(3x3)ありますが,回転の自由度は3なので,回転行列以外の方法でパラメータ化しなければいけません.今回は,Angle-Axisを用いてパラメータ化しています.
daily-tech.hatenablog.com
X.Sympy を使った微分計算
ヤコビアン,ヘッシアンともに近似をつかわず(!)計算しています(笑)さすがに手計算はむりだったので,Sympyを使っています.
daily-tech.hatenablog.com
X.勾配ベクトルとヘッシアンの作成
バンドル調整における典型的なヘッシアンの構造を踏まえて,勾配ベクトル/ヘッシアンを作ります.
daily-tech.hatenablog.com
X.直線探索
最急降下法・ニュートン法を実行する上で必要となる直線探索のアルゴリズムを調べています.
daily-tech.hatenablog.com
X.ニュートン法を用いた解法
力ずくで計算したヘッシアンを使って,最適化をやってみました.ただ,簡単には収束せず...
X.レベンバーグ・マルカート法を用いた解法
非線形最小二乗法のデファクトスタンダードである,レベンバーグマルカート法を試しました.
X.Ceres-Solver を用いた解法
最後に,,,同じ問題設定を Ceres-Solver を用いて解きました.ヤコビアンもヘッシアンも計算せず,あっという間ですね...
X.参考文献
Bundle Adjustment – A Modern Synthesis
SBA: A software package for generic sparse bundle adjustment
Numerical Optimization (Springer Series in Operations Research and Financial Engineering)
- 作者:Nocedal, Jorge,Wright, Stephen
- 発売日: 2006/06/01
- メディア: ハードカバー
- 作者:健一, 金谷
- 発売日: 2005/09/01
- メディア: 単行本
X.実験に使ったソースコード
ソースコードはまだちょっと整理中でして,個々のエントリのアップする時には間に合わせるようにします.
github.com
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
をモットーに,いろんな実験をしていけたらと思います!ではでは,今年もよろしくお願いします!
ビックマック
前を向きたいとき,僕はビックマックを食べる.
仕事が終わらず,2.5㎡のオフィスでコードとにらめっこしていた時,僕はビックマックを食べた.
クリスマスに誰も出迎えてくれなかったときも,僕はビックマックを食べた.
「年末,休めないんじゃね..」と不安になりながらも,僕はビックマックを食べた.
この先も,何度もくじけそうになるだろう.だから,これからもよろしく.
ということで,仕事終わったーーー!!!なんとか冬休みを迎えられそうです.
仕事納めのご褒美として,久々にビックマックを食べてみました.「ご褒美にマックって,どんな貧乏社長やねん」って言われそうですが,お金的にというよりは,カロリー&背徳感的なご褒美ですかね.
結局今年はあんまりブログ書けませんでしたが,,,,冬休みはちょっと頑張ってアップできればと思います(゚Д゚)ノ
オフィス
どうもどうも,社長です(笑).
近況
ちょっと更新が滞ってしまったんですが,実は7月頭から会社が始動しておりまして,せっせと働き始めてます.だいぶスタートダッシュをかまさないといけない状況でありまして(笑),ブログの更新頻度がちょっと落ちてしまうと思うのですが,ちょこちょことアップしますのでまた暇なときに見ていただければありがたいですm(__)m
(今年はつくばチャレンジも不参加なので,SFMもなんとかやり切ります!)
オフィス契約
勢いで法人を立ち上げ,一人会社なので自宅を登記したのですが,これ,国税庁のホームページで住所出てくるんですよね...登記後にそれに気づいたんですが,後の祭りでして...会社のホームページ立ち上げたり,ブログで会社名を公表しようと思ってたんですが,身に覚えのないピザとかがいっぱい届いたら困るので(笑)オフィスを借りて,住所を登記しなおすことにしました.
で,まともなオフィスはさすがに高いので,レンタルオフィスを借りることにしました.残念ながら南流山近辺にはなかったので,ビズサークルという会社が運営している市川妙典のレンタルオフィスにしました.南流山からは少し離れているので,マンション更新のタイミングで市川市に引っ越そうと思います.ということで,なかなか港区民への道のりは遠そうです(笑)
肝心のオフィスですが,共益費込みで月額 32000 円でして,何とかやりくりできそうな額です.下にオフィスの写真を載せてます.自分でいうのもなんですが,あまりの ”おままごと” っぷりに自然と笑みがこぼれます(爆)が,今をときめく Apple もガレージから始まったことを考えると,先はどうなるかわかりませんよ~.目指せ100兆円起業 ( ゚Д゚)
訪問調査
ある程度の規模の会社になると,初回取引では相手法人の信用調査をします.ただ,法人調査部門を会社毎に持つことは難しいので,普通は調査会社に外注します.日本の場合だと,東京商工リサーチと帝国データバンクの2社が有名どころだと思います.
実は,会社(=自宅)に東京商工リサーチの調査員の方が来ました!訪問調査というよりは,家庭訪問みたいな形になってしまったんですが,一応となりのコンビニでおーいお茶を購入して,丁重にもてなしました.ということで僕の会社,企業情報データベースに乗ってます!スコアはおそらく最低レベルだと思いますが(笑)
それでは,今後ともわが社をよろしくお願いしますm(__)m
(まだ先になりそうですが,会社住所の登記変更が完了したら,社名&ホームページも公開します!)
Fisher Vector を用いた類似画像検索
ということで,ちょっと時間がかかってますが類似画像検索まで来ました!
X. 類似画像検索
類似画像検索とは何か?という話なんですが,Googleの画像検索のように,クエリ画像に対して「よく似ている」画像を検索することです.
「類似画像検索がなぜSFMに必要なのか?」という話なんですが,SFMを実行する場合,ペアとなる写真には同じ物体が写っている必要があります.例えば東京タワーを3次元復元したい場合,東京タワーが写っていない写真を使っても意味がありません.もし手持ちの写真に東京タワーが写っていない写真も含まれているとしたら,これを除く必要があります.下記の例では,スカイツリーの写真を除外しています.
復元したい物体のスケールが「東京タワー」とか「奈良の大仏」程度なら,頑張れば一枚の写真に全てを収めることができます.が,「ローマの町並み」とか「国際通り」とかになってくると,写真一枚に全てを収めることは当然できないので,オーバーラップのある写真をペアとして繋いでいく必要があります.
上の写真は国際通りのお菓子御殿(笑)ですが,3つの写真を使ってお菓子御殿を復元したい場合,左の写真は少なくとも中央の写真と,右の写真も少なくとも中央の写真とペアにならなければ3枚の写真間につながりができません.ということで,類似している画像を検索することによって,クエリ画像とペアになる写真を見つけていきます.(※写真が時系列に並んでいる場合,おそらく直近の数枚の写真がペアになるので,類似検索は必ずしも必要無いです.)
X. Fisher Vector とは
いつもの如く TheiaSfM で使われていた方法を参考に, Fisher Vector による方法を用いました.Fisher Vector の原理を理解するには,統計とか情報幾何学とかをかじる必要があって,数学力の低い私は今頭を抱えているのですが...このへんもある程度理解できたらまたエントリを残せればと思います.
原理は難しいものの,ライブラリを使えば処理フローはそれほど複雑ではなかったので,簡単に書き留めておきます.大まかに4つの処理からなっていまして,
1.全画像の全SIFT特徴量を用いて,GMM(混合ガウスモデル)を生成.
ー>生成過程では,全画像のSIFT特徴量を入力してやればOKです.多すぎる場合,間引いた方がいいかもしれません.この混合ガウスモデルが特徴量の生成モデルになります.
2.生成した GMM を用いて,個々の画像を Fisher Vector に変換.
ー>一つ一つの画像データを Fisher Vector に変換します.ここではFisher Vector が一つの画像を表現しています.今回の設定では一つの画像から大体15000個位の SIFT 特徴量が生成されたのですが,各識別子の次元が 128 なので,SIFT を用いて一つの画像を表現すると,15000 x 128 = 1920000 位の次元のベクトルになります.が,これでは大きいので,生成した GMM のパラメータ(の勾配)空間を使って次元削減してやります.この次元削減によって,今回のケースだと一枚の画像が 8000 次元くらいのベクトル(= Fisher Vector)で表現できるようになりました.
3.全画像ペア間の Fisher Vector の距離を計算します.
ー>「2つの画像が似ている」=「2つの画像の Fisher Vector が似ている」ということになります.あとは,全ペアの Fisher Vector の距離を総当り計算して行列として保持しておけば,画像間の類似関係をいつでも参照することができます.
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) をダウンロードして使いました.
クエリ画像としては,下記の画像を使いました.
で結果ですが,スコアの低い(類似度が高い)画像から順に上位10枚の画像が下記のようになりました.下記の結果を見てもわかるのですが,やっぱり天候とか昼夜の条件によってスコアに開きが出てしまうのは仕方ないのかもしれませんね...この辺の影響を抑える方法があれば,つくばチャレンジにも使えるようになるかも.
スコアの分布は下記のようになります.最低スコア(類似度が高い)0.516687を皮切りに,上限2.3くらいまで連続的に分布していきます.
X. 実験に使ったコード
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 の計算のチュートリアルは下記になります.
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. 実験に使ったコード
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
特徴点(SIFT)のシリアライズとその他もろもろ
画像間類似度を計算する方法を調査してるんですが,画像の枚数が多くなってくると毎回SIFTを計算するのも辛くなってくるので,結果を保存する方法を調べることにしました.このトピックに限らずですが,参考にしたのはTheiaSfMというオープンソースプロジェクトです.githubのリンクは最後に載せてます.
X.Cereal
boost とか,いくつかシリアライズのためのライブラリがあるみたいなんですが TheiaSfM で Cereal が使われていたこともあって,ひとまず使ってみることにしました.ドキュメントも用意してくれていて,とても親切な感じですm(_ _)m
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 と STL の Vector を使って表現します.
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
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; }