最近傍探索（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
中部大 藤吉先生の説明スライドですが，めっちゃわかりやすいです．
Scale Space
ということで，下記２つのエントリを経てやっとこさ Scale Space のとこまで来ました．
daily-tech.hatenablog.com
daily-tech.hatenablog.com
で， 微分フィルタ及び二階微分フィルタ（Laplacian Filter）のエントリの最後で，
「スケールが 2σ 位の被写体が写っている画像に対して，標準偏差が σ のガウスフィルタと二階微分フィルタを適用すると，被写体の部分が強く反応する．」
という話をしました．このエントリではこれを実際に実験＆確認して，最後に Scale Space の定義を紹介します．
Laplacian of Gaussian Filter を適用してみる．
ここでは，実験用に下記の画像を使います．
この画像に写っているひまわりですが，大体下記のようなサイズ感になっています．
ここで，赤の円が半径 45pix，青の円が半径 10pix，緑の円が半径 5pix になります．ちなみに，ひまわりは花の部分が黄色，種の部分が黒なので，円で囲ってある部分のエッジはx方向からスキャンしていくと，
「下がりエッジー＞上がりエッジ」
になります．なので，二階微分は
「下がりエッジー＞上がりエッジと上がりエッジー＞下がりエッジ」
の組み合わせになるので，「強く反応する＝円の中心付近で出力値が大きくなる」ことを意味します．
それではこの写真に対して，前回のエントリで実験した 「Laplacian of Gaussian Filter」を適用して結果を見てみます．理論的にはガウシアンフィルタの標準偏差が10pix くらいになったとこで青の円が強く反応（＝明るくなって）して，45pix くらいになったとこで赤の円が強く反応（＝明るくなって）するはずですが，，，実験してみます．
σ = 1.6
σ = 2.2
σ = 3.2
σ = 4.5
σ = 6.4
σ = 9.0（青の円付近のひまわりが強く反応．）
σ = 12.8
σ = 18.1
σ = 25.6（赤の円付近のひまわりが強く反応．）
σ = 36.2（赤の円付近のひまわりが強く反応．）
σ = 51.2
σ = 72.4
ちょっと写真でいっぱいになってしまいましたが，σが大きくなるにつれて強調されるひまわりが画像遠方から近方に写ってくる様子がわかると思います．遠方はほぼ黄色（花びらの部分）と茶色（種の部分）で構成されているのに対し，近方になってくると緑色（葉や茎）の要素が出てくるため，これによってエッジ構造が複雑になり，強め合い・弱め合いが起こることによって実際のサイズとピークを取るσからの乖離が大きくなってしまいました．が，ある程度は理論に沿った動きをすることがわかりました．
Scale Spaceの定義．
上記ひまわり畑の実験を踏まえると，LoG Filter (Laplacian of Gaussian Filter) のσを変化させることで，スケールが大体 σ 位の物体・構造を検知することができるようになりました．そこで，画像の縦・横に加え，新たにσの次元を加えることによって，画像中にどんなスケールの物体が写っているかがある程度わかるようになります．下記，簡単ですがいつものポンチ絵です．
図中にスケールがそれぞれ r1, r2, r3 である３つの物体が書かれていますが，σ を変化させて LoG filter をかけ，その反応が強くなった箇所を書き留めておくことによって，書き留めた座標 (x, y) を中心とする物体のスケールが大体 σ くらいであるということがわかるようになります．SIFTではこの原理を使って特徴点を抽出します．ということで，次は SIFT に入ります．あー．時間がかかる．．．
実験に使ったコード．
import os import cv2 import numpy as np def scale_space_test1(): cur_dir = os.path.dirname(os.path.abspath(__file__)) path = cur_dir + '/../data/himawari.jpg' img = cv2.imread(path) gray_img = cv2.cvtColor(img, cv2.COLOR_RGB2GRAY) k = np.sqrt(2) sigma = 1.6 # for sigma in sigmas: for i in range(12): sigma_3_radius = 10 * sigma * 2 # Applying Gaussian filter gauss_kernel_1d = cv2.getGaussianKernel( ksize=int(sigma_3_radius), sigma=sigma) gauss_kernel_2d = np.outer( gauss_kernel_1d, gauss_kernel_1d.transpose()) log_filter = cv2.Laplacian(gauss_kernel_2d, cv2.CV_32FC1, ksize=31) laplacian_of_gaussian = cv2.filter2D( gray_img, cv2.CV_32FC1, log_filter) # Normalize value for visualization. min_val, max_val, _, _ = cv2.minMaxLoc(laplacian_of_gaussian) abs_max = max(abs(min_val), abs(max_val)) tmp = (laplacian_of_gaussian * (0.5 / (2.0 * abs_max)) + 0.5) * 255.0 laplacian_of_gaussian_uchar = np.uint8(tmp) title = 'Sigma = ' + str(sigma) cv2.namedWindow(title) cv2.imshow(title, laplacian_of_gaussian_uchar) #cv2.imwrite(title + str('.png'), laplacian_of_gaussian_uchar) cv2.waitKey(0) cv2.destroyWindow(title) sigma = k * sigma if __name__ == "__main__": scale_space_test1()
参考文献．
Scale-space theory: A basic tool for analyzing structures at different scales
Distinctive image features from scale-invariant keypoints
微分フィルタ及び二階微分フィルタ（Laplacian Filter）
前回の Gaussian Filter のエントリに引き続き，微分フィルタについてもエントリを残しておこうと思います．前回のエントリは下記の通り．
daily-tech.hatenablog.com
一階微分フィルタ・二階微分フィルタとは？
その名の通り，隣接する画素値の変化率を計算し，その値を出力とするフィルタです．一般的に画像のエッジを取りたいときによく使われます．こちらもググればいっぱい有益な情報が出てくるので，ここでは実際に使ってみてどんな出力が得られるか見てみたいと思います．
２D画像に一階微分フィルタ・二階微分フィルタを適用する．
手始めに，下記のサンプル画像に対して一階微分フィルタ・二階微分フィルタを適用します．下記の画像はわが母校．．．ではありませんが，東工大の大岡山キャンパスです．この画像に対して，x方向の一階微分フィルタ，ｙ方向の一階微分フィルタ，二階微分フィルタを適用してみます．
原画像
ｘ方向の一階微分フィルタを適用
ｘ方向に関して，上がりエッジがある場所では値が高く（白く），下がりエッジがある場所では値が低く（黒く）なっている様子がわかります．
ｙ方向の一階微分フィルタを適用
ｙ方向に関して，上がりエッジがある場所では値が高く（白く），下がりエッジがある場所では値が低く（黒く）なっている様子がわかります．
二階微分（ラプラシアン）フィルタを適用
二階微分フィルタの適用結果ですが，一階微分フィルタの応答とは異なり，よく見るとエッジがあるところが白黒になっていることがわかると思います．
二階微分（ラプラシアン）フィルタの結果拡大図
上の画像の一部拡大版ですが，二階微分（ラプラシアン）フィルタを適用すると，エッジの存在する箇所が「黒白，もしくは白黒」になっていることがわかります．
１D画像に微分フィルタを適用する．
ここではそれぞれのフィルタの動きを見るために，問題を簡単化して１D画像で考えます．前回のエントリ同様に白帯の画像を考えます．上から順に，「原画」「原画のｘ方向断面」「ｘ方向の一階微分フィルタ出力」「二階微分フィルタの出力」となっています．
上記の図から，二階微分フィルタの出力がなぜエッジ付近で「黒白，もしくは白黒」になることがわかるかと思います．上がりエッジが立つと，そのエッジを達成するために一階微分（勾配）が
「０ー＞正の値ー＞０」
と変化し，この勾配を達成するために二階微分が
「０ー＞正の値ー＞負の値ー＞０」
となるためです．二階微分フィルタを適用した大岡山キャンパス画像の出力結果ともつながりました．
ひとまず，一階微分フィルタ・二階微分フィルタの概要はこのくらいにしておきます．ここではエッジがとてもくっきりしていて，ある意味ではステップ入力のような白帯を使いました．次に，この白帯をGaussian Filter でぼかしたときの二階微分（ラプラシアン）フィルタの出力変化を見てみます．
Gaussian Filter と二階微分（ラプラシアン）フィルタの関係
下記，原画に対して Gaussian Filter をかけたときの一階微分フィルタ，二階微分フィルタの応答です．σが大きくなるに従って勾配が緩やかになり，一階微分フィルタ，二階微分フィルタの応答も緩やかになっている様子がわかると思います．
σの標準偏差を持つガウスフィルタを適用した場合，ステップ入力の勾配が小さく（緩やかに）なるため，一階微分・二階微分の値は小さくなっていきます．この出力だと形がわからないので，一階微分フィルタ・二階微分フィルタの出力を正規化して，どんな形になっているのか見てみます．（一階微分の出力値にσを，二階微分の出力値にσ^2をそれぞれかけてあげると正規化できます．）
正規化した画像から，おもしろいことがわかります．４段目の二階微分フィルタの値を見てほしいのですが，白帯の上がりエッジ・下がりエッジにて発生した波がガウスフィルタの標準偏差（σ）を大きくしていくに従って互いに接近しはじめ．．．．σ = 16 で完全に合体しました．ちなみに，この白帯の幅は 32pix なのですが，白帯の幅のちょうど半分くらいの標準偏差を持つガウスフィルタをかけてあげると，二階微分フィルタの出力が極値を迎えることがわかりました．ということで，
「スケールが 2σ 位の被写体が写っている画像に対して，標準偏差が σ のガウスフィルタと二階微分フィルタを適用すると，被写体の部分が強く反応する．」
ことがわかります．今回は一つの画像に対して σ を変化させて微分フィルタの応答を見ましたが，この σ の方向の次元を 「Scale Space」 と呼びます．Scale Space の極値を見ることによって，（＝今回実験したように，σを変えながら二次微分フィルタの極値を探ることによって）被写体の大体のスケールがわかることになります．SIFTではキーポイント検出でこれを使っているのですが，それはまた次のエントリで．．．．
実験に使ったコード．
参考文献．
Scale-space theory: A basic tool for analyzing structures at different scales
Distinctive image features from scale-invariant keypoints
Gaussian Filter を用いた画像の暈し．
「SFMをやります！」と言いながらだいぶ基礎的なエントリになるのですが，
「SFMをやる！」
ー＞「特徴点はSFMの基礎だから，SIFTとAKAZEくらいは抑えておこう．」
ー＞「まてよ，Scale Space ってなんや？」
という流れで，フィルタリングを簡単にまとめておくことにしました．「SFMシリーズ！」としてエントリをまとめてみようと思ったのですが，そうするとだいぶ先まで見通してエントリを書かないといけないので，小出しに個別のエントリを書くことにしました．一通りやりきることができたら，あとでリライトします．
ということで，国際通りの３次元復元まで，道のりは遠いですね．．．
Gaussian Filter とは．
テンプレートマッチングや機械学習するときに，画像の前処理としてフィルタリングをすることもあると思います．自分はこのあたりあまり考えず，「ちょっと画像暈したいなあ．」とかって具合に使ってたんですが（笑），Gaussian Filter には意味があったみたいです．Gaussian Filter はググればいっぱい有益な情報が出てくるので，かっちりとした話はそちらを見ていただくとして，ここではその意味合い的な部分を掘り下げられればと思います．
参考サイト１
imagingsolution.blog.fc2.com
標準偏差 (σ)と暈しの程度を考える．
Gaussian Filterはガウスカーネルを持つフィルタなので，フィルタの設定値として標準偏差（σ）を指定します．OpenCV の関数でいうと，，，
cv2.getGaussianKernel(ksize, sigma, ktype)
の sigma の部分です．この関数を使う上で気を付けないといけないことは，sigma を大きくした場合，それに伴って ksize の値も大きくしないといけないということです．（σが大きくなるにつれて，当該ピクセルの値を計算するために必要な隣接ピクセルの数が増えていくため．）
それでは，画像に Gaussian Filter を適用する場合にこの σ がどんな意味を持つのかちょっと考えてみます．実験として，下記画像に対してガウシアンフィルタを適用してみます．まず簡単のために，x方向だけの1Dガウシアンフィルタを使いたいので，画像はとてもシンプルにy方向には何も変化がないものにしています．
ちなみに，上記の実験画像は高さが61pix, 幅が501pixで，画像の中心の白帯の幅が11pix になります．次に，σを変化させたときの1Dガウシアンカーネルをプロットしてみると，下記のようになります．
で，実際に実験画像に上記の1Dガウシアンカーネルを適用した結果が下記になります．ここでは「画像の中で幅11pixの領域として映っているものが，σ を変化させたときにどのようにぼかされていくか？」を見ています．左側のグラフは，右側画像をx軸に平行な断面で切って横から見た図です．
上図を見るとわかると思うのですが，σ が11を超えたあたりから急速に白帯が消失してく様子がわかります．σが 11 を超えるまでは，どちらかというとエッジの部分がぼかされていたのに対して，11を超えてくると黒画像領域に白帯が溶け出してきます．ということでガウスフィルタを使うときは，
「画像のなかで暈したいスケール (pix) を決め，σをそのスケールよりも少し大きくとる」
ことをすれば，いい具合に画像をぼかせそうです．
画像にGaussian Filterを適用する．
次に，実際の2D画像に Gaussian Filter を適用してみたいと思います．ここでは実際に上の議論が成立しているかどうかを見るために，下記のような実験画像を用意しました．四角形，丸，三角形が並んでいますが，上から順に 10 x 10pix，20 x 20pix, 40 x 40pix, 80 x 80pix, 80 x 160pix の大きさで書いてます．一番下はちょっと扁平にしてみました．
で，この画像に対して，「それぞれの図形スケール + α」の標準偏差を持つガウスフィルタを適用してみます．ここまでの議論が正しければ，図形のスケールよりも大きい標準偏差を持つガウスフィルタを適用すると図形は消えるはず！といことで，やってみます．
σ = 5
１番スケールの小さい図形（10 x 10 pix）はだいぶ暈されてますが，まだ残ってます．
σ = 180
５番目スケールの図形（80 x 160 pix）が消えました．
ということで，思ったような結果になりました．ただ，暈しなので完全には消えておらず薄っすらと残ってますが，大体のスケールを決めてやれば暈しの効果をある程度操ることはできそうです．ということで，次は Laplacian FIlter と Scale Space の話をします．
実験に使ったコード．
参考文献．
Scale-space theory: A basic tool for analyzing structures at different scales
Distinctive image features from scale-invariant keypoints
春休み！
2020年もはや5月になってしまいましたが，，，推定5万人の僕のブログ読者の皆様（笑），元気にされてますでしょうか！！ちょっとブログの更新タイミングを逃してしまい，だいぶ間が空いてしまったんですが，令和元年初のポストです．
１．ここまでの報告
昨年８月から始まったフリーランス生活，ここまでいい感じに来ています．といってもまだ９か月ですが，，，ちょうど先月いっぱいで１件目のお客さんの仕事を離れたのですが，現場の雰囲気が良かったこともあり，とても充実した９か月でした．
今までのキャリア的に「R＆D畑で論文を読み漁りながらスクラッチで物をつくりまくる．」というよりは，どちらかというと「量産開発やある程度形になっている物を改編・導入」という仕事が多かったので，アルゴリズム検討＆スクラッチ開発してアウトプットが出せるかどうかビクビクでした．が，数年にわたってロボット作ったり画像処理の勉強をしてきた複利効果が表れてきたのか，一応それなりにアウトプットを出すことができ，お客さんにも喜んでもらえたと思っています．
まだまだ修行の身ではありますが，テーマを絞って継続的に勉強なり実験なりすれば，２～３年くらいたてば実感できる程度にはスキルが上がってくるということに気づけたこともとても大きな収穫でした．
２．やりたいこと！
実はこれからちょっと春休みします（笑）
本当は「ハワイ or 八重山諸島」だったのですが，コロナ的にそれはとても難しそうなので，南流山で勉強＆モノづくりに勤しみます．３か月くらい春休みしようと思っているんですが，想定以上に景気が怪しくなってきたので，，，様子を見ながらちょっと短くするかもです．
で，肝心のやりたいこと！なんですが，下記になります．いくつかは去年のブログにも書いていたと思うのですが，脳力不足のせいでどれも未完になってます．．．今年はどれか一つでも終わればいいな．．．
１．つくばチャレンジ完走！
ことしで６回目の出場になる，「つくばチャレンジ」の完走です．が，，，コロナのせいで今年は実験走行・本走行が中止となってしまい．．．今年は完走できたはずなので（笑）とても残念ですが，2021年に備えます(｀・ω・´)ゞ
子分との勉強会もちょっと要調整ですね．
２．国際通りの三次元復元．
那覇の国際通りをSFMで三次元復元したい！
３．ディープラーニングの勉強．
自分のスキルセットに「機械学習・深層学習」がまだない状態でして．．．学習系のスキルの有り無しで受けられる仕事の数がだいぶ変わってくるので，そろそろ勉強したい．
４．数学の勉強．
ロボティクスにしても画像処理にしても最近の論文はかなり最適化を使っていて，数学のバックグラウンドがないとあんまり理解が深まらないので，すこし遠回りかと思いつつも数学（主に線形代数・最適化）の勉強をしたいと思っています．
５．会計・金融の勉強．
２０１９年に人生で初めて青色申告したんですが，いままで工学系のことしか勉強してこなかった自分にとってはこれがとても新鮮で．
３．２０２０年，やること！
ということでどれも捨てがたいのですが，さすがに３ヵ月しかないので浅く広くといっても絞らざるを得ず．．．ひとまず，下記の３つにします．
・国際通りの三次元復元．
３度目の正直です（笑）SFMやります！ただ，今すぐ沖縄には行けないので，南流山駅前でやってみて，うまくいく目途が立ったらコロナの収束具合を見て那覇にいきます．こいつは引き続きオープンソースとにらめっこして，適宜下記を読みながらやっていきます．
Multiple View Geometry
Multiple View Geometry in Computer Vision
- 作者:Hartley, Richard,Zisserman, Andrew
- 発売日: 2004/03/25
- メディア: ペーパーバック
・数学の勉強．
実は，数学に関しては年始から少しずつやり始めているんですが，下記を課題図書にしてます．もともとそんなに数学が強くないことと，中身が濃すぎてまとめ記事を書くのがすごく難しそうなんですが，このあたりってそれほどブログ記事になってるものが少ないような気がするので，できればエッセンスを抜き出して記事にしたいと思っています．（これは今年中目標です（笑））
Matrix Computation
Matrix Computations (Johns Hopkins Studies in the Mathematical Sciences)
- 作者:Golub, Gene H.,Van Loan, Charles F.
- 発売日: 2012/12/27
- メディア: ハードカバー
Convex Optimization
Convex Optimization (English Edition)
- 作者:Boyd, Stephen,Vandenberghe, Lieven
- 発売日: 2004/03/08
- メディア: Kindle版
Numerical Optimization
Numerical Optimization (Springer Series in Operations Research and Financial Engineering)
- 作者:Nocedal, Jorge,Wright, Stephen
- 発売日: 2006/06/01
- メディア: ハードカバー
ちなみに，なんで洋書やねん ( ﾟДﾟ) って話があると思うんですが，アマゾンのレビューがすこぶるいいことと，全部ネットに pdf で落ちてるからです！
・会計・金融の勉強．
で，もう一つ！独立っていうとなんだか大げさですが，一応青色申告をする身分になったので，会計の基礎的なことを勉強したいと思います．ただ，勉強するだけだと面白くないので，企業の決算報告書とかをある程度読めるようになります！で，真剣に株式投資をやってみます！（＜－これが真のモチベーションなのかも（笑））ということで，下記，課題図書です．
- 作者:國貞克則
- 発売日: 2016/10/13
- メディア: 新書
- 作者:國貞克則
- 発売日: 2016/10/13
- メディア: 新書
財務3表実践活用法 会計でビジネスの全体像をつかむ (朝日新書)
- 作者:國貞克則
- 発売日: 2012/07/13
- メディア: 新書
![会社四季報ワイド版 2020年2集春号 [雑誌]](https://m.media-amazon.com/images/I/61h1F0LxWoL._SL160_.jpg "会社四季報ワイド版 2020年2集春号 [雑誌]")
- 発売日: 2020/03/16
- メディア: 雑誌
投資で一番大切な20の教え 賢い投資家になるための隠れた常識
- 作者:ハワード・マークス
- 発売日: 2012/10/23
- メディア: 単行本
- 作者:リンチ,ピーター,ロスチャイルド,ジョン
- 発売日: 2001/03/01
- メディア: 単行本
こちらの課題図書はすべてお金を出して購入させていただきました（笑）
あとは，日経新聞のオンラインと四季報オンラインを申し込むかどうか．．．．迷い中です．
ということで，２０２０年はここから怒涛のブログ更新率を見せつけますので，５万人の読者のみなさま，どうぞよろしくおねがいしますm(__)m
