オフィス
どうもどうも，社長です（笑）．
近況
ちょっと更新が滞ってしまったんですが，実は７月頭から会社が始動しておりまして，せっせと働き始めてます．だいぶスタートダッシュをかまさないといけない状況でありまして（笑），ブログの更新頻度がちょっと落ちてしまうと思うのですが，ちょこちょことアップしますのでまた暇なときに見ていただければありがたいですm(__)m
（今年はつくばチャレンジも不参加なので，SFMもなんとかやり切ります！）
オフィス契約
勢いで法人を立ち上げ，一人会社なので自宅を登記したのですが，これ，国税庁のホームページで住所出てくるんですよね．．．登記後にそれに気づいたんですが，後の祭りでして．．．会社のホームページ立ち上げたり，ブログで会社名を公表しようと思ってたんですが，身に覚えのないピザとかがいっぱい届いたら困るので（笑）オフィスを借りて，住所を登記しなおすことにしました．
で，まともなオフィスはさすがに高いので，レンタルオフィスを借りることにしました．残念ながら南流山近辺にはなかったので，ビズサークルという会社が運営している市川妙典のレンタルオフィスにしました．南流山からは少し離れているので，マンション更新のタイミングで市川市に引っ越そうと思います．ということで，なかなか港区民への道のりは遠そうです（笑）
肝心のオフィスですが，共益費込みで月額 32000 円でして，何とかやりくりできそうな額です．下にオフィスの写真を載せてます．自分でいうのもなんですが，あまりの ”おままごと” っぷりに自然と笑みがこぼれます（爆）が，今をときめく Apple もガレージから始まったことを考えると，先はどうなるかわかりませんよ～．目指せ１００兆円起業 ( ﾟДﾟ)
訪問調査
ある程度の規模の会社になると，初回取引では相手法人の信用調査をします．ただ，法人調査部門を会社毎に持つことは難しいので，普通は調査会社に外注します．日本の場合だと，東京商工リサーチと帝国データバンクの２社が有名どころだと思います．
実は，会社（＝自宅）に東京商工リサーチの調査員の方が来ました！訪問調査というよりは，家庭訪問みたいな形になってしまったんですが，一応となりのコンビニでおーいお茶を購入して，丁重にもてなしました．ということで僕の会社，企業情報データベースに乗ってます！スコアはおそらく最低レベルだと思いますが（笑）
それでは，今後ともわが社をよろしくお願いしますm(__)m
（まだ先になりそうですが，会社住所の登記変更が完了したら，社名＆ホームページも公開します！）
Fisher Vector を用いた類似画像検索
ということで，ちょっと時間がかかってますが類似画像検索まで来ました！
X. 類似画像検索
類似画像検索とは何か？という話なんですが，Googleの画像検索のように，クエリ画像に対して「よく似ている」画像を検索することです．
「類似画像検索がなぜSFMに必要なのか？」という話なんですが，SFMを実行する場合，ペアとなる写真には同じ物体が写っている必要があります．例えば東京タワーを３次元復元したい場合，東京タワーが写っていない写真を使っても意味がありません．もし手持ちの写真に東京タワーが写っていない写真も含まれているとしたら，これを除く必要があります．下記の例では，スカイツリーの写真を除外しています．
復元したい物体のスケールが「東京タワー」とか「奈良の大仏」程度なら，頑張れば一枚の写真に全てを収めることができます．が，「ローマの町並み」とか「国際通り」とかになってくると，写真一枚に全てを収めることは当然できないので，オーバーラップのある写真をペアとして繋いでいく必要があります．
上の写真は国際通りのお菓子御殿（笑）ですが，３つの写真を使ってお菓子御殿を復元したい場合，左の写真は少なくとも中央の写真と，右の写真も少なくとも中央の写真とペアにならなければ３枚の写真間につながりができません．ということで，類似している画像を検索することによって，クエリ画像とペアになる写真を見つけていきます．（※写真が時系列に並んでいる場合，おそらく直近の数枚の写真がペアになるので，類似検索は必ずしも必要無いです．）
X. Fisher Vector とは
いつもの如く TheiaSfM で使われていた方法を参考に， Fisher Vector による方法を用いました．Fisher Vector の原理を理解するには，統計とか情報幾何学とかをかじる必要があって，数学力の低い私は今頭を抱えているのですが．．．このへんもある程度理解できたらまたエントリを残せればと思います．
原理は難しいものの，ライブラリを使えば処理フローはそれほど複雑ではなかったので，簡単に書き留めておきます．大まかに４つの処理からなっていまして，
１．全画像の全SIFT特徴量を用いて，GMM（混合ガウスモデル）を生成．
ー＞生成過程では，全画像のSIFT特徴量を入力してやればOKです．多すぎる場合，間引いた方がいいかもしれません．この混合ガウスモデルが特徴量の生成モデルになります．
２．生成した GMM を用いて，個々の画像を Fisher Vector に変換．
ー＞一つ一つの画像データを Fisher Vector に変換します．ここではFisher Vector が一つの画像を表現しています．今回の設定では一つの画像から大体15000個位の SIFT 特徴量が生成されたのですが，各識別子の次元が 128 なので，SIFT を用いて一つの画像を表現すると，15000 x 128 = 1920000 位の次元のベクトルになります．が，これでは大きいので，生成した GMM のパラメータ（の勾配）空間を使って次元削減してやります．この次元削減によって，今回のケースだと一枚の画像が 8000 次元くらいのベクトル（= Fisher Vector）で表現できるようになりました．
３．全画像ペア間の Fisher Vector の距離を計算します．
ー＞「２つの画像が似ている」＝「２つの画像の Fisher Vector が似ている」ということになります．あとは，全ペアの Fisher Vector の距離を総当り計算して行列として保持しておけば，画像間の類似関係をいつでも参照することができます．
４．クエリ画像とその他画像の距離チェックし，距離が近いものを類似画像とする．
ー＞ クエリ画像のインデックスをもつ 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) をダウンロードして使いました．
クエリ画像としては，下記の画像を使いました．
で結果ですが，スコアの低い（類似度が高い）画像から順に上位１０枚の画像が下記のようになりました．下記の結果を見てもわかるのですが，やっぱり天候とか昼夜の条件によってスコアに開きが出てしまうのは仕方ないのかもしれませんね．．．この辺の影響を抑える方法があれば，つくばチャレンジにも使えるようになるかも．
スコアの分布は下記のようになります．最低スコア（類似度が高い）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
で，今回は実行ファイルが３つありまして，順を追って実行する必要があります．
１．作成済のSIFT特徴量を用いて，GMMの生成
（flagfile.txt に SIFT 特徴量，画像ファイルのディレクトリパスを指定し，生成した GMM モデルのパラメータファイルの出力パスを指定します．）
./build/train_gmm_model --flagfile=./flagfiles/train_gmm_model.txt
２．生成したGMMを用いて，全画像の Fisher Vector を計算．各画像間の距離マトリクスを作成する．
（flagfile.txt に SIFT 特徴量のディレクトリパス，生成した GMM のパラメータファイルパスを指定し．生成するマッチングマトリクスの出力ファイルパスを指定します．）
./build/create_matching_matrix --flagfile=./flagfiles/train_gmm_model.txt
３．クエリ画像のファイル名を指定すると，類似している上位画像を出力します．
（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; }
最近傍探索（Nearest Neighbor Search）
ということで，なんだかんだ会社設立やその他もろもろが思った以上に大変で，，，ちょっと止まってたんですが，ひとまず小休止に入ったので引き続き SFM をやりたいと思います！
X. 最近傍探索とは．
まず，最近傍探索とは何ぞや？って話なんですが，Wikipedia の定義を見ると，
「距離空間における最も近い点を探す最適化問題の一種」
とあります．「距離空間」には「マンハッタン距離」とか「ユークリッド距離」とかいろいろとありますが，ここでは簡単のために「ユークリッド距離」前提としてエントリを進めていきます．で，これだけ見ても意味わからないので，まず簡単に一次元の例を考えてみたいと思います．
Brute Force Macher １次元データの最近傍探索
この問題は，「クエリ点（=4.2）にもっとも近い値を探索対象の中から探す．」ということになりますが，上記は力任せ方による解法です．このやり方だとクエリ点と探索対象点のすべての距離を計算して，距離が最も近いものを選択するという流れになります．このやり方では，比較対象点がn個存在する場合，探索回数が O(n) になります．なので，比較対象点が１億点あるときは1億回の比較が必要になります．
Binary Treeによる１次元データの最近傍探索
次に，Binary Tree を用いた最近傍探索になります．力任せ探索とは異なり，比較対象点が n 個存在する場合，探索回数が O(log n) のオーダになります．実際には数値のばらつき具合によりますが，うまく行けば比較対象点が１億点あるときは，約27回程度の比較になります．めちゃ減りましたね．上図では，ツリーの葉ノードが探索対象を表しており，四角のノードが分岐のための閾値を表しています．閾値は単純に当該ノード以下の探索対象の平均値を用いています．
で，今回のエントリでは，「点が１次元じゃなくて，高次元のユークリッド空間にある場合はどうやって探索すんの？」という話になります．高次元になったときでも力任せ探索は有効で，クエリ点に対して一番近い点を確実に見つけてくれます．ただ，時間がかかります．一方で， Binary Tree を高次元に適用しようとすると，ちょっと工夫が必要そうです．
X. SFMの処理フローと計算コスト．
で，過去に書いた SFM エントリのちょっとしたおさらいになるんですが，細かな処理を除くと SFM の大きな処理フローとは下記になります．
１．各写真から特徴点の抽出．
２．視点が似ている（＝同じようなものが被写体になっている）写真のペアを見つける．
【探索コスト大】全ての写真ペアを比較し，視点が似ているかチェック．
３．見つけた写真のペアに紐づく特徴点をマッチング．
【探索コスト大】ペアの全ての特徴点を比較して，最近傍ペアを推定する．
４．カメラの位置と点群を推定する．
で，２と３のところで探索の話が出てきます．２では写真ペア，３では特徴点ペアということで，ちょっと探索の内容は異なりますが，SFMの復元対象物体が大きくなるほどこの２と３にかかるコストが大きくなっていきます．奈良の大仏とかその位のスケールの構造物でも写真の枚数がある程度大きくなるはずなので，２と３を工夫しないと計算量が爆発します．今回のエントリでは，３の特徴点探索の話をします．
※ちなみに那覇の国際通りは1.6キロあるので，２と３を工夫するだけだとまだ不十分で，Hybrid SfM か，写真の時系列順序を使う方法を試したいと思ってます．
X. SIFTを用いた場合の SFM の特徴点マッチング．
今回は特徴点に SIFT を使おうと思っていますが， SIFT では画像中で見つかったキーポイントに一つづつに識別子（Descriptor）が付随しています．
つまり，生成したキーポイントの比較に識別子を使うことになるのですが，この識別子は上図にあるように128個のビンを持つヒストグラムになっています．つまり，「128次元のユークリッド空間内に存在する点の最近傍探索」になります．
X. 今回解剖した３つの手法．
一つ一つのエントリがそこそこの分量になりそうだったので，別のポストにします！書き終わったら順次リンクを貼っていきます．
Brute Force Matching
エントリを書き終わったらリンクしますm(_ _;)m
KD Tree Search
エントリを書き終わったらリンクしますm(_ _;)m
Cascade Hashing for 3D Search
エントリを書き終わったらリンクしますm(_ _;)m
X. ３つの手法による最近傍探索の計算時間．
計算時間を全く考慮に入れなければ Brute Force Matching で良いのですが，そういうわけにもいかないので，計算時間とマッチング精度のトレードオフを取ることになります．で，今回は実験を簡単にするために，同じ画像から抽出したキーポイントと識別子を比較しました．ただ，比較対象となる特徴点セットの中に必ず完全一致する特徴点が見つかります．これでは実際に SFM において使う用途と状況が異なってしまうので，探索対象データに少しずつノイズを混ぜていって，探索時間や精度がどう変化するか見てみました．
実験に用いた画像
取得できた特徴点数：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. 実験に使ったコード
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.
代表取締役
ということで，なっちゃいました（笑）
ブログの更新が滞ってしまったんですが，実は絶賛会社設立中です．もともと法人化に興味はあったのですが，「個人事業がうまく行って，いつか会社を設立できればいいなあ．．．」くらいに考えていたところ，ちょうど仕事の話をいただきまして，法人格でないと発注できないということだったので「良い機会かも！」と思ってトライしてみることにしました．
ただ，会社設立までの具体的な作業や書類等，何一つ把握してなかったので，そこから作業をはじめました．
情報調査のために購入した書籍
個人事業主からの法人化になるのですが，下記の書籍を購入しました．一通り必要なことをまとめてくれていて，分かりやすかったです．
改訂4版 個人事業を会社にするメリット・デメリットがぜんぶわかる本
- 作者:俊輔, 関根
- 発売日: 2019/08/13
- メディア: 単行本
- 作者:寺内正樹
- 発売日: 2010/10/19
- メディア: 単行本
会社設立に使ったサービス
上記の本をちゃんと読んだあとに下記のサービスを使って法人化のための書類一式を作成したのですが，必要事項を入力していけば，ほぼ修正無しで会社設立できました．．．．ただ，事前に一連の流れを知っておくほうが良いとは思いますが．．．
www.freee.co.jp
直近のドタバタ．．．
・会社設立を決意！Amazonで本を買う．
・Amazon本で一通りの流れを調査．
・TODOリストを起こす．（法人登記までの流れ，銀行口座作成 etc etc...）
・市役所で個人印鑑の登録＆印鑑証明取得．
・freee の会社設立サービスを使って，書類一式の作成．
・会社名考える．．．
・会社名決定，定款最終版作成．
・定款認証＠公証役場，会社の印鑑発注．
・印鑑が届き，法務局で登記申請．
・登記完了！法人印鑑証明，登記簿謄本を取得．
・銀行への説明資料作成．．．
・税務署へ法人設立届出書を出した後，銀行へ口座開設申請＆資料送付．
ということで，頂いた仕事の話が急ぎだったこともありだいぶ突貫工事だったんですが，なんとか会社設立までたどりつけそうです．ネットで調べたところ法人用の銀行口座開設は結構難しいらしいのですが，，，頑張って資料作成して泣き落としできればと（笑）
飛躍の５月
なんだかこのところ，５月は毎年いろいろな転機が訪れる月になっています．
一昨年の５月は，「ドイツ人と日本人のはざまで苦悩すべく，ドイツへ２ヶ月の旅．」
daily-tech.hatenablog.com
去年の５月は，「個人事業主になることを決意．退職し，沖縄へ（笑）」
daily-tech.hatenablog.com
今年の５月は，「代表取締役として，会社を取り締まる立場に（笑）」
今後
ちょっと酔っていることもあって，だいぶ軽いノリのエントリになってしまいました．．．エンジニアとしても（当然のことながら会社としても）まだまだ未熟者ですが，「コンピュータビジョンとロボティクス」をメインテーマとして地道にやっていきたいと思います．会社名とか会社のホームページとか，もうちょっと形になってきたらブログでも紹介できればと思います．それでは，今後共よろしくお願いします！
SIFT ~ 実装編（全体構成）
SIFTの理論と実装のエントリ，どうやってまとめようか悩んでたんですが，ちょっと絵を書き始めると時間がかかりすぎて前に進めなくなってしまうので，ひとまず大雑把にまとめておくことにしました．コードはだいぶリファクタして読みやすく（自分としては．．．）なったはずなので，コードを追いかけてもらうとSIFTの実装がどんな感じになっているかはざっくりと解ると思います．一通りSFMに目処がついた段階で，各エントリをリライトできればと思います．
あと，ブログの最後に参考文献のリンクを貼っておいたのですが，中部大の藤吉先生が作られているスライドがとっても分かりやすかったので，かっちりとした理論の部分はそちらを参考にしていただければと思います．
X．全体構成 特徴点抽出
OpenCVの cv::xfeatures2d::SIFT クラスをそのままパクってリファクタしたのですが，論文の構成にできるだけ近づくようにリファクタしました．論文の流れと一致するよう，関数を「特徴点検出」と「記述子作成」の２つに分けたのですが，ここでは特徴点検出の部分を見ています．
１．CreateInitialImage
処理の前段階として，
① 入力画像を Float 型 グレー画像へ変換．
② 入力サイズの２倍に画像を拡大．この画像が First Octave (first_octave_idx =-1) となります．
③ ２で拡大した画像に対して，ガウシアンフィルタを適用します．
を実施します．
２．Octave 数の計算．
つぎに，必要となる Octave 数を計算します．処理画像（①で前処理された画像）に対して，高さと幅の大きい方を 2 で割っていったときに，一辺の長さが 4pix を下回らないところまで画像を縮小していきますが，この値が Octave 数になります．この値は３で計算される Gaussian Pyramid の段数と同じになります．
３．BuildGaussianPyramid
画像ピラミッドを作成します．ここでは処理画像をピラミッドの底辺として，Octave 数段のピラミッドを作成します．
４．BuildDifferenceOfGaussianPyramid
３で生成した Gaussian Pyramid の各 Octave の画像間で差分を計算し，Difference of Gaussian Pyramid を作成します．
６．RemoveDuplicateSorted
重複しているキーポイントを削除します．
７．RetainBest
特徴点の検出上限数が指定されている場合，スコアが低いものを取り除いて上限数に収まるようにします．
８．CompensateInitialScalingToOriginalImage
１のステップで入力画像サイズを２倍にされているので，計算されたキーポイントのスケール，位置，サイズを求まった値に対して半分にします．
該当部コード
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 ループの入れ子構造になってますが，本質的な計算をしているのは下記の部分になります．
１．キーポイント周辺領域計算．
特徴点抽出の部分で抽出されたキーポイントに対して，記述子を計算するための周辺領域サイズを計算します．
２．必要なバッファ準備．
ヒストグラム計算に必要なバッファを確保します．必要なものとしては，「Mag : 勾配強度」「Ori : 勾配角度」「W : 重み（中心点からの距離による）」「RBin : 最終記述子の列方向成分」「CBin : 最終記述子の行方向成分」「hist : ヒストグラム計算領域」の６つです．
３．ComputeGradientForPrePatch
キーポイント周辺領域を分割した小領域に対して，勾配，勾配方向を計算します．
４．UpdateHistogramBasedOnTriLinearInterpolation
３では小領域に対する勾配，勾配方向を計算しましたが，これを記述子に統合します．記述子は１６領域（＝４ｘ４）からなり，一つ一つの領域に８ビンのヒストグラムをもちます．小領域の計算結果がそれぞれこの１６ヒストグラムの各ビンに振り分けられていきます．
５．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
中部大 藤吉先生の説明スライドですが，めっちゃわかりやすいです．
