2018年11月30日金曜日

Colmapによるカメラの内部パラメータ推定

カメラキャリブレーションをして内部パラメータが欲しいときは, OpenCVのcameraCalibration というサンプルプログラムはあまりいい選択肢ではない. マーカの検出が不安定でキャリブレーション用のパターン全体が画像内に写っていなければ, データを取得することができず画像の隅の方までマーカを写すことができないから.

Colmapという三次元復元 (SfM + MVS) ソフトは,カメラ位置や特徴点の三次元位置だけなく内部パラメータまで推定してくれるので, カメラキャリブレーションに使用することができると聞いたので早速使ってみた. 似たようなソフトにBundlerがあるが,複数の画像を入力したとき,複数の画像で内部パラメータが独立に推定されるのと, カメラモデルを選択することができずカメラキャリブレーションには向いていない. Colmapでは,OpenCVのカメラモデルやその他のモデルを選択でき,特徴点により計算できるので, キャリブレーション専用のパターンやボードを準備する必要がない上, 画像の隅々まで

入力画像と歪み補正結果

今回のキャリブレーションで左が入力画像の1枚,右がそのレンズ歪み補正結果.直線であるはずの歪んだ構造が歪み補正結果で直線に戻っていたら成功. そういう意味ではなかなかうまく補正されているように見える.

本来の使い道は三次元復元なのでその結果ももちろん得られる.

実行手順

  1. 動画(video.mp4)から連番画像(img-xxx.png)ファイルを作る
    ffmpeg -i video.mp4 -f image2 -vcodec png img-%03d.png
    
  2. COLMAP-dev-windows-legacy-cuda.zipダウンロード+解答
  3. Colmap.bat起動
  4. File → New project
    1. Newクリックして作業用フォルダ選択→適当な名前.dbを入力して保存
    2. Selectクリックして画像のフォルダ選択
    3. Saveクリック
  5. Processing → Future extraction
    1. Colmap (https://colmap.github.io/)でcamera modelをFULL_OPENCVに設定
    2. Camera modelのオプションShared for all images (全ての画像が同じ内部パラメータ)
    3. Camera modelのオプションCustom parameters(初期値)に640,640,640,512,0,0,0,0,0,0,0,0を設定
    4. Extractクリック
  6. Processing → Future matching → Exhaustive (すべての画像のペアでマッチング)
    1. Runクリック
  7. Reconstruction → Reconstruction options → Bundle → Camera parameters
    1. refine_***を全部チェック
  8. Reconstruction → Start reconstructionクリック (約1時間)
  9. Reconstruction → Bundle adjustment
    1. refine_****を全部チェック
    2. Runクリック
  10. File → Export model as text(ここでキャリブレーションデータ抽出終わり)
  11. Reconstruction → Dense reconstruction(ここからはおまけ.メッシュモデルができる)
    1. Undistortionクリック
    2. Stereoクリック
    3. Fusionクリック
    4. Poissonクリック
    5. Delaunayクリック

パラメータ初期値の決め方

今回,「Custom parameters(初期値)に640,640,640,512,0,0,0,0,0,0,0,0を設定」したけど,これはどうやって与えたのか? COLMAPのカメラモデルのうち「Full OpenCV camera model」を用いたので,「$f_x, f_y, c_x, c_y, k_1, k_2, p_1, p_2, k_3, k_4, k_5, k_6$」の順に値が並んでいる. これらは,OpenCVのパラメータなので「Camera calibration With OpenCV」を見て欲しい.

今回は,入力画像の水平方向の解像度が1280画素で,カメラの画角がおおよそ90度くらいだったので, 「画角と解像度から内部パラメータ行列を作る」を参考にすると, camera matrixの部分は決まる.レンズ歪みのパラメータについては全部ゼロ(つまり歪み無し)を初期値として, Colmapの最適化の過程で最適な値を求めてもらった.

Colmapを閉じた後の再ロード方法

  1. Colmapを開く
  2. File→Open project
  3. プロジェクトファイルを開く(Export model as textで出力されたもの)
  4. File→Import model
  5. 「Export model as text」で出力したフォルダを選択

参考

2018年11月22日木曜日

標準偏差と標準誤差を直感的に理解する

この記事は,研究室のゼミでせっかくグラフをエラーバー付きで示してプレゼンしたのに, もしくは104.3±4.5のように誤差つきで計測値を論文に書いたのに, それって標準偏差それとも標準誤差?などと聞かれた人向け(それ以外にも誤差の示し方は半値幅とか2σとかいくらでもある). 標準偏差の意味は,なんとなくデータのばらつき具合で, 精度(もしくは誤差)を表すのに都合が良い,というところまでは理解しているけど, よく似た名前でよく似た数式の標準誤差とどうちがうのか,いまいち教科書を読んでも良くわからない,という人向け. 「三角形は,コンビニのおにぎりの形」くらいのざっくりした説明をする.

標準偏差と標準誤差の違い

まず名前が違う

当たり前だけど名前が違う. 名前が物事を分かりにくくするということが多々あって, 私もこの名前のせいで理解しようとする気持ちが長年失せ続けた人間の一人. まず,共通部分の「標準」だけど,はっきり言って「標準」ということに殆ど意味はないので無視した方がいいい. 「統一規格の」というくらいの意味なのだろうか?命名された経緯は良くわからない.

次に「偏差」と「誤差」という部分が違いだけど, これも2種類の統計量の違いを的確にあらわしている訳ではないので,無視して丸暗記するしかない. しいて自分で名前をつけるなら,標準偏差は「ばらつき幅」,標準誤差は「平均値精度」が良いと思っている. そう,わざわざ全然違う名前にするということは,意味が違うということ. この記事では名前の分かりにくさが理解を妨げてるという立場をとってるので,「標準偏差(ばらつき幅)」のように常にセットで用いる.

数式が違う

それぞれの定義式を説明しはじめると,ややこしくなるので両者の関係式から. 教科書を読むとごちゃごちゃ数式がでてきて嫌になってしまって, どうやって計算するのかが書かれている最後の部分だけ見るとこんな感じになっている.

\[ 標準誤差(平均値精度)=\frac{標準偏差(ばらつき幅)}{\sqrt{データ数}} \]

要するに標準偏差が計算できれば標準誤差は簡単に計算できるということ. 「というか,ほとんど同じやん」と思ったのでは? 全然違う名前をつけた割りには,標準誤差(平均値精度)は単に標準偏差(ばらつき幅)をデータ数のルートで割っただけのもの. 要するに標準偏差(ばらつき幅)が大きなときは標準誤差(平均値精度)も大きくなる,という意味では同じようなもの.

こう考えると全然意味が分からなくなる. この結果は,もともと異なる統計量を計算した偶然の結果であってこの関係詞期から両統計量の違いを理解しようとすると理解できなくなる. 実際教科書では,標準偏差は直感的なイメージとともになぜそういう定義がされたか説明されるが,標準誤差は中心極限定理を持ち出して説明される.

使い方が違う

標準偏差
データのばらつきを表すのに用いる.データ単体を示して云々語るときはこちら.
標準誤差
データから得られた平均値の精度を表すのに用いる.2つ以上のデータセットの平均値の差を比較する場合はこちら.

以上です.詳しくは教科書を読んで勉強して欲しい.

参考

2018年10月4日木曜日

シンクライアント上でUnity Proの認証を自動化する

目的

大学の演習室などに導入されるシンクライアント上でUnity Proアカウントを認証させ,Unityを起動後に特別な操作を必要とせず使える状態にしたい.

環境・条件・制約

大学の演習室と言っても色々あると思う.自分が経験した環境は次のようなものだった.

  • Windows 10
  • Unity 2017.3.1f1 (64bit)
  • OSイメージがネットワーク経由で提供されるタイプのシンクライアント.
  • イメージの書き換えは容易ではない.
  • 教員が提供する共有フォルダを学生アカウントから読み取りのみ可能.

方法概略

Unity Proアカウントを認証させて起動するとログイン情報とライセンスファイルが生成されることが, 新規作成されるファイルを丹念に調べると分かった. ログイン情報は他のPCでもそのままコピーして使用できるけど, ライセンスファイルにはPCの情報がハッシュ化されて保存されるので 他のPCのファイルを使いまわすことができない.そこで,次の方法で必要なファイルを収集して, シンクライアントをログインするたびに共有フォルダなどから必要なファイルをコピーしてからUnityを起動する.

  • 各クライアントで下の方法で「ログイン情報を収集」し,ファイルを配置してログインする.
  • シリアルキーを入力してUnity Proアカウントを認証させた上で,「ライセンスファイルを収集」する.
  • シンクライアント上でUnityを実行するときはバッチファイルなどで「ログイン情報とライセンスファイルを配置」した上でUnityを起動する.

ログイン情報を収集

%username%はアカウント名に置き換わる.このログイン情報は一定期間で期限切れになるみたいなので,定期的に更新する必要がある.

copy /Y C:\Users\%username%\AppData\LocalLow\Unity\Browser\Cookies d:\

ライセンスファイルを収集

%COMPUTERNAME%はマシン名に置き換わる. これで異なるクライアントのライセンスファイルをUSBメモリやネットワークドライブの同じフォルダに集める. このファイルはどうも期限がなさそう. PCの情報が含まれているので演習を受ける人が勝手に持ち返っても利用できない.

copy /Y C:\ProgramData\Unity\Unity_lic.ulf D:\%COMPUTERNAME%\Unity_lic.ulf

ログイン情報とライセンスファイルを配置

この内容をバッチファイルに入れて,Unity使用時に実行するようにすればいい.

copy /Y d:\Cookies C:\Users\%username%\AppData\LocalLow\Unity\Browser\Cookies

copy /Y d:\Cookies C:\Users\%username%\AppData\LocalLow\Unity\Browser\Cookies
C:\ProgramData\Unity\Unity_lic.ulf

2018年8月23日木曜日

畳み込みのフーリエ変換に出現する$2\pi$

2つの関数$g(x), h(x)$の畳み込み積分$g*h$のフーリエ変換$\mathcal{F}[g*h]$の公式をみると色々とバリエーションがあって混乱するので整理してみた.

同じ畳み込みなのに公式いろいろ

\[ \def\diff{{\rm d}} \] \begin{eqnarray} g * h&=&\int^{\infty}_{-\infty} g(y) \, h(x-y) \diff y \end{eqnarray}

畳み込み積分がこのように定義されているとすると,そのフーリエ変換は次のように,各関数のフーリエ変換の単純な積になる. ただし,資料によっては積の前に$2\pi$がついていたり,ついていなかったり. 明らかにこれらは異なるものなので,ちょっと混乱してしまう.

\begin{eqnarray} \mathcal{F}[g * h]&=&\ &\mathcal{F}[g] \, \mathcal{F}[h]\\ \mathcal{F}[g * h]&=&\ \sqrt{2\pi} &\mathcal{F}[g] \, \mathcal{F}[h]\\ \mathcal{F}[g * h]&=&\ 2\pi&\mathcal{F}[g] \, \mathcal{F}[h] \end{eqnarray}

結論だけいうと,この違いはフーリエ変換の定義の違いからくるもの. 後で示してるようにこの違いの原因には逆変換は関与しなくて順変換だけが影響してくる.つまり, 両辺のフーリエ変換の順変換の回数が異なる(左辺1回,右辺2回)ため,フーリエ変換の順変換側についている係数の積が等しくなるように$\sqrt{2\pi}$などの係数が出現しているだけと考えたら覚えやすい.

フーリエ変換の定義いろいろ

そのフーリエ変換の定義はなんと3つもある(もっとあるかも).

\begin{eqnarray} \mathcal{F}[g] &=& &\int^{\infty}_{-\infty} g(x) e^{-i\omega x} \diff x,\ \ & g(x) &=&\frac{1}{2\pi}&\int^{\infty}_{-\infty} \mathcal{F}[g](x) e^{i\omega x} \diff \omega\label{def1}\\ \mathcal{F}[g] &=& \frac{1}{\sqrt{2\pi}}&\int^{\infty}_{-\infty} g(x) e^{-i\omega x} \diff x,\ \ & g(x) &=&\frac{1}{\sqrt{2\pi}}&\int^{\infty}_{-\infty} \mathcal{F}[g](x) e^{i\omega x} \diff \omega\label{def2}\\ \mathcal{F}[g] &=& \frac{1}{2\pi}&\int^{\infty}_{-\infty} g(x) e^{-i\omega x} \diff x,\ \ & g(x) &=& &\int^{\infty}_{-\infty} \mathcal{F}[g](x) e^{i\omega x} \diff \omega\label{def3}\\ \end{eqnarray}

なぜこんな異なる定義があるのかというと,フーリエ変換した後,逆変換して,もとに戻ってくれば問題ないから. それさえ満たせば係数はどのようにでも定義可能で,どれが正しいというものでもない. だから,フーリエ変換にまつわる公式を使う場合は,自分がどの定義を使っているのか認識して,一貫して同じ定義から導出される公式を使わないとおかしなことになる.工学系の教科書では一番上の定義が多い,らしいけど,どれも結構な頻度で現れるし,ネット上の資料だとどの定義から出てきたものか書いてないことが多い.

周波数としての代わりに波数$f=\omega/2\pi$が用いられる場合

これもややこしいんだけど,周波数空間を表す変数として,位相角ではなくて,$2\pi$分の波の数で表す場合がある. この場合もフーリエ変換の定義式が異なり次のようになる.

\begin{eqnarray} \mathcal{F}[g] &=& &\int^{\infty}_{-\infty} g(x) e^{-i2\pi f x} \diff x,\ \ & g(x) &=&&\int^{\infty}_{-\infty} \mathcal{F}[g](f) e^{i2\pi f x} \diff f\\ \mathcal{F}[g] &=& \frac{1}{\sqrt{2\pi}}&\int^{\infty}_{-\infty} g(x) e^{-i\omega x} \diff x,\ \ & g(x) &=&\frac{1}{\sqrt{2\pi}}&\int^{\infty}_{-\infty} \mathcal{F}[g](\omega) e^{i\omega x} \diff \omega\\ \mathcal{F}[g] &=& &\int^{\infty}_{-\infty} g(x) e^{-i\omega x} \diff x,\ \ & g(x) &=&\frac{1}{2\pi}&\int^{\infty}_{-\infty} \mathcal{F}[g](\omega) e^{i\omega x} \diff \omega\\ \end{eqnarray}

ただし,この場合は水平軸の単位が異なるだけなので,垂直軸方向の振幅自体に影響はなくて,上とまったく同じ使い分けをすれば良い.

どこでそんな違いが出てくるのか,証明までさかのぼってみた.

$2\pi$が出現する瞬間

フーリエ順変換の前についている係数が$1$だったり$\frac{1}{\sqrt{2\pi}}$だったり$\frac{1}{2\pi}$だったりして全部確認するのは面倒なので,これらの係数をまとめて全部${\color{red}A}$とおいてみる.$\def\A{{\color{red}A}}$

\begin{eqnarray} \A = 1,\ \frac{1}{\sqrt{2\pi}},\ {\rm or} \frac{1}{2\pi} \end{eqnarray}

目立つように大文字かつ赤字で.フーリエ順変換の定義が登場するたびにこの${\color{red}A}$がでてくるのでそこに注意.

\begin{eqnarray} %\def\diff{{\rm d}} \mathcal{F}[g * h] \ &=&\ \A \int^{\infty}_{-\infty} \left( \int^{\infty}_{-\infty} g(y) \, h(x-y) \diff y \right) \, e^{-i\omega x} \diff x \\ &=&\ \A \int^{\infty}_{-\infty} g(y) \left( \int^{\infty}_{-\infty} h(x-y) \, e^{-i\omega x} \diff x \right) \diff y \\ &=&\ \A \int^{\infty}_{-\infty} g(y) \left( \int^{\infty}_{-\infty} h(t) \, e^{-i\omega(t+y)} \diff t \right) \diff y \\ &=&\ \A \int^{\infty}_{-\infty} g(y) \left( \int^{\infty}_{-\infty} h(t) \, e^{-i\omega t} \diff t \right) \, e^{-i\omega y} \diff y \end{eqnarray}

外側のフーリエ順変換の中に更にフーリエ順変換っぽい形が見えてきたけど係数$\A$が無いのでここではまだフーリエ順変換の定義とはみなさないでおく.ここで両辺に$\A$を掛ける.

\begin{eqnarray} \A~\mathcal{F}[g * h]&=&\ \A \int^{\infty}_{-\infty} g(y) \left( \A \int^{\infty}_{-\infty} h(t) \, e^{-i\omega t} \diff t \right) \, e^{-i\omega y} \diff y \\ &=&\ \A \int^{\infty}_{-\infty} g(y) \mathcal{F}[h] \, e^{-i\omega y} \diff y \\ &=&\ \A \int^{\infty}_{-\infty} g(y) \, e^{-i\omega y} \diff y \ \mathcal{F}[h]\\ \mathcal{F}[g * h]&=&\ \A~\mathcal{F}[g] \, \mathcal{F}[h] \end{eqnarray}

はい,ついに出てきました!この通り順変換の前の係数だけに依存して畳み込みのフーリエ変換の公式が変化することが分かりました.めでたしめでたし.

2018年8月19日日曜日

フーリエ変換によるODOGフィルタ計算の高速化

Oriented Difference of Gaussians (ODOG)という視覚モデル[1]を用いた画像のフィルタリングをフーリエ変換を用いて高速化する方法を考えてみる. ODOGは,明度知覚に関する様々な効果をうまく説明できるモデルとして知られていて, ODOGモデルは1枚の画像が与えっれたとき,人間が知覚する明度が画素値になった1枚の画像に変換する. この変換の一部は畳み込み演算になっているので,高速フーリ変換アルゴリズムを使うと 単純に畳み込み演算で計算するよりも,画像やODOGフィルタのカーネルを一旦フーリ変換し, 結果を逆フーリエ変換して元に戻す方がむしろ高速になる.

ODOGって

まず,ODOGがどういう数式モデルなのかを説明しないと,話が始まらない. 簡単に言うと,このODOGは,DoG[2]という異なるサイズのガウス関数の差で,二次微分のような効果があるフィルタの二次元版のフィルタがベースになっている. 二次元版のDoGは,方向に依存したフィルタになっていて,様々なスケールと方向のDoGを,元の画像にかけて,ある重みで足し合わせたものになる.

具体的には,入力画像$p(x_1,x_2)$に対して,ある方向$\alpha=k\pi/3$ $(k=0,\cdots,5)$, あるスケール$\sigma_1=3.0(1/2)^{l}$ $(l=0,\cdots,7)$のDoGカーネル$f_{k,l}(x_1, x_2)$を畳み込んだものが$q_{k,l}(x_1, x_2)$だとすると, ODOGの出力$R(x_1,x_2)$は次の式で定義される.

\begin{eqnarray} \label{eq:1} R(x_1, x_2)=\frac{1}{6}\sum_{k=0}^5 \frac{Q_k(x_1, x_2)}{||Q_k(x_1, x_2)||_2} \end{eqnarray} \begin{eqnarray} \label{eq:2} Q_k(x_1, x_2)=\sum_{l=0}^6 w_l\cdot q_{k,l}(x_1, x_2) \end{eqnarray} \begin{eqnarray} w_l = (8/3\sigma_1)^{-\frac{1}{10}} \end{eqnarray} \begin{eqnarray} \label{eq:4} q_{k,l}(x_1, x_2)=f_{k,l}(x_1, x_2)*p(x_1,x_2) \end{eqnarray}

ただし,$||\cdot||_2$は,$L^2$ノルムであり,次のように定義される.

\begin{eqnarray} \left( ||Q_k(x_1, x_2)||_2 \right)^2 = \int_\Omega Q(x_1, x_2)^2 dx_1 dx_2 \end{eqnarray}

この式のもう少し詳しい説明は,論文[1]か[3]を参照して欲しい.

畳み込み演算をフーリエ変換で高速化する

上の式(\ref{eq:1})から(\ref{eq:4})の中でどの部分をフーリエ変換で高速化できるのか? フーリエ変換で畳み込み演算が掛け算になることを利用できるのは式(\ref{eq:4}). そこで,まず式(\ref{eq:4})をフーリエ変換する.

関数$f$のフーリエ変換を$\mathcal{F}[f]$と表すと,式(\ref{eq:4})両辺のフーリエ変換は次のようになる.

\begin{eqnarray} \label{eq:6} \mathcal{F}[q_{k,l}]=\mathcal{F}[f_{k,l}]\ \mathcal{F}[p] \end{eqnarray}

ODOGのカーネル$f$のフーリエ変換$\mathcal{F}[f_{k,l}]$は, 「ODOGカーネルのフーリエ変換」 で紹介しているように,予め解析的に計算することができる.

足し算を分解する

フーリエ変換は線形写像なので足し算の部分は分解でき,係数は前に出すことができる. 具体的には,式(\ref{eq:2})の総和記号以下は足し算なので次のように計算できる.

\begin{eqnarray} \mathcal{F}[Q_k]&=&\mathcal{F}\left[\sum_{l=0}^6 w_l\cdot q_{k,l} \right]\\ &=&\sum_{l=0}^6 w_l\cdot \mathcal{F}\left[q_{k,l} \right] \end{eqnarray}

ここで式(\ref{eq:6})を代入すると,変数$l$を含まない$\mathcal{F}[p]$は係数と同じように扱えるので総和記号の前に出すことができる.

\begin{eqnarray} \mathcal{F}[Q_k]&=&\sum_{l=0}^6 w_l\cdot \mathcal{F}[f_{k,l}]\ \mathcal{F}[p]\\ &=& \mathcal{F}[p] \sum_{l=0}^6 w_l\cdot \mathcal{F}[f_{k,l}]\\ \end{eqnarray}

この$\sum_{l=0}^6 w_l\cdot \mathcal{F}[f_{k,l}]$の部分は,$k$を定めると予め数値計算しておくことができる. 変数$k$の値は0から7の8つのみなので,予め計算したものを配列として保存しておけば実行時計算する必要はない.

また,$\mathcal{F}[p]$については,$p(x_1,x_2)$が画像なので,この画像に対してフーリエ変換を一度計算すれば良い.

正則化項の計算を高速化する

$\mathcal{F}[Q]$の計算法は示したけど,正則化項$||Q_k(x_1, x_2)||_2$の計算をするには 単純には一旦$\mathcal{F}[Q]$を逆フーリエ変換$\mathcal{F^{-1}}$して$Q$を求める必要がある. しかし,パーシバルの定理によるとフーリエ変換前後のエネルギーは(係数を除いて)保存されるので,実は逆フーリエ変換は必要ない.

\begin{eqnarray} E_k=\int_\Omega Q_k(x_1, x_2)^2 dx_1 dx_2 = \frac{1}{2\pi}\int_\Omega \mathcal{F}[Q_k](t_1, t_2)^2 dt_1 dt_2 \end{eqnarray}

再び足し算を分解する

式(\ref{eq:1})の両辺をフーリエ変換して,足し算部分を分解し,上で計算したものを代入する.

\begin{eqnarray} \label{eq:12} \mathcal{F}[R]=\frac{1}{6}\sum_{k=0}^5 \frac{\mathcal{F}[Q_k]}{\sqrt{E_k}} \end{eqnarray}

高速化アルゴリズム

フーリエ変換を高速フーリエ変換(FFT)に変えると,一連の処理は次のようになる.

  1. 入力画像$f$に対してFFTすることにより$P=\mathcal{F}[p]$を計算する.
  2. $F_k = \sum_{l=0}^6 w_l\cdot \mathcal{F}[f_{k,l}]$を計算して配列に蓄積する.
  3. $\mathcal{F}[Q_k]$を計算する.
  4. $\mathcal{F}[Q_k]$の全画素の自乗和により$E_k$を計算する.
  5. $\mathcal{F}[R]$を計算する.
  6. 逆高速フーリエ変換により$R$を計算する.

参考文献

  1. Barbara Blakeslee and Mark E. McCourt: "A multiscale spatial filtering account of the White effect, simultaneous brightness contrast and grating induction," Vision Research, vol. 39, issue 26, pp. 4361-4377, 1999.
    doi: 10.1016/S0042-6989(99)00119-4
  2. "Difference of Gaussians" in Wikipedia
    https://en.wikipedia.org/wiki/Difference_of_Gaussians
  3. Barbara Blakeslee, Davis Cope, and Mark E. McCourt1: "The Oriented Difference of Gaussians (ODOG) Model of Brightness Perception: Overview and Executable," Behavior Research Methods, pp. 306–312, vol. 48, no. 1, 2016.
    doi: 10.3758/s13428-015-0573-4

2018年8月17日金曜日

ODOGカーネルのフーリエ変換

ここではOriented Difference of Gaussians (ODOG)という視覚モデル[1]のフーリエ変換を手計算して確認してみる. ODOGは,明度知覚に関する様々な効果をうまく説明できるモデルとして広く知られている. ODOGモデルは1枚の画像が与えっれたとき,人間が知覚する明度を画素値として持つ画像に変換しフィルタとして機能する数式モデル. このモデルは,初等関数の簡単な組み合わせでできているので,簡単にフーリエ変換できる. フーリエ変換するメリットは,このフィルタを入力画像に畳み込む計算を,フーリエ変換すると単なる掛け算になって,高速化できること. とりあえずそうした背景は一切忘れて,ODOGのカーネルをフーリエ変換してみる.

ODOGのカーネル関数って?

その肝心のODOGのモデルだけど,全ての数式を説明するっだけでそれなりのスペースが必要になるので,ここではODOGのカーネル関数だけをフーリエ変換してみる. ODOGのカーネル関数には,角度とスケールのパラメータがあり,様々な角度とスケールのカーネル関数でフィルタした結果を足し合わせて最終的な画像を出力する仕組みになっている. そのカーネル関数$f(y_1, y_2)$を数式で表すと次のようになる.

\[ \def\infint{\int_{-\infty}^\infty} \] \begin{eqnarray} \label{eq:1} f(y_1, y_2) & = & \frac{1}{2\pi\sigma_1^2}e^{-\frac{y_1^2 + y_2^2}{2\sigma_1^2}} - \frac{1}{2\pi\sigma_1\sigma_2}e^{-\frac{1}{2}(\frac{y_1^2}{\sigma_2^2} + \frac{y_2^2}{\sigma_1^2})}\\ % %\label{eq:2} %&=&\frac{1}{2\pi\sigma_1}e^{-\frac{y_2^2}{2\sigma_1^2}} %\left(\frac{1}{\sigma_1}e^{-\frac{y_1^2}{2\sigma_1^2}} - \frac{1}{\sigma_2}e^{-\frac{y_1^2}{2\sigma_2^2}}\right)\\ \end{eqnarray}

ただし,画素の位置を表す変数は$(x,y)$や$(u,v)$ではなくて,論文[2]に合わせて$(y_1,y_2)$にした.

この式は,異なる2つのガウス関数の差に見える. 実際に$y_1$-軸方向の断面を見るために$y_2=0$を代入すると次のように2つの異なるガウス関数の差になっている.

\begin{eqnarray} \label{eq:2} f(y_1, 0) & = & \frac{1}{2\pi\sigma_1^2}e^{-\frac{y_1^2}{2\sigma_1^2}} - \frac{\sigma_2}{\sigma_1}\frac{1}{2\pi\sigma_2^2}e^{-\frac{y_1^2}{2\sigma_2^2}} \end{eqnarray}

一方,$y_2$-軸方向の断面を見るために$y_1=0$を代入すると次のように1つの異なるガウス関数になっている.

\begin{eqnarray} \label{eq:3} f(0, y_2) & = & \frac{1}{2\pi\sigma_1^2}e^{-\frac{y_2^2}{2\sigma_1^2}} - \frac{1}{2\pi\sigma_1\sigma_2}e^{-\frac{y_2^2}{2\sigma_1^2}} = \left(1 - \frac{\sigma_1}{\sigma_2} \right) \frac{1}{2\pi\sigma_1^2}e^{-\frac{y_2^2}{2\sigma_1^2}} \end{eqnarray}

このように,普通のDoG[3]と異なり,方向をもつ構造になっている.

画像に対する二次元フーリエ変換

教科書でよく出てくるのは一次元の信号$f(t)$なんかのフーリエ変換. $f$のフーリエ変換を$\mathcal{F}[f]$と表すと,周波数$\omega$を用いて次のように定義されている.

\[ \mathcal{F}[f](\omega)=\infint f(t)\ e^{-i \omega t}d\omega \]

画像に対する二次元フーリエ変換を計算するには,$x$方向と$y$方向(正確にはここでは$y_1$方向と$y_2$方向)の一次元のフーリエ変換を二重にかければよい. フーリエ変換は積分なので二重積分になる. 周波数空間($s_1,s_2$)でのカーネル関数を$F(s_1,s_2)$とすると次のようになる.

\begin{eqnarray} \label{eq:4} F(s_1, s_2) = \infint \left(\infint f(y_1, y_2) e^{-i y_1 s_1}dy_1\right) e^{-i y_2 s_2}dy_2 \end{eqnarray}

ここでも論文[2]に合わせて周波数の変数は$(s_1,s_2)$とした.

足し算のフーリエ変換は分けることができる

式(\ref{eq:4})の()の中の一次元フーリエ変換をまず計算してみる. そのために式(\ref{eq:1})を$f(y_1, y_2)$に代入する. 式(\ref{eq:1})は,$a-b$の形になっている.フーリエ変換の大事な性質は線形性. つまり$a$と$b$それぞれのフーリエ変換を$A$と$B$だとすると,$a-b$のフーリエ変換は単純に$A-B$になるので, 式(\ref{eq:1})のフーリエ変換は次のように各項単独のフーリエ変換の差になる.

\begin{eqnarray} \label{eq:6} \infint f(y_1, y_2) e^{-i y_1 s_1}dy_1 &=& \infint \frac{1}{2\pi\sigma_1^2}e^{-\frac{y_1^2 + y_2^2}{2\sigma_1^2}} e^{-i y_1 s_1}dy_1\\ &-& \infint \frac{1}{2\pi\sigma_1\sigma_2}e^{-\frac{1}{2}(\frac{y_1^2}{\sigma_2^2} + \frac{y_2^2}{\sigma_1^2})}\ e^{-i y_1 s_1}dy_1\\ %------------------------------------------------------------------------------- &=& \label{eq:8} \frac{1}{2\pi\sigma_1^2}e^{-\frac{y_2^2}{2\sigma_1^2}} \infint e^{-\frac{y_1^2}{2\sigma_1^2}}\ e^{-i y_1 s_1}dy_1\\ &-& \label{eq:9} \frac{1}{2\pi\sigma_1\sigma_2} e^{-\frac{y_2^2}{2\sigma_1^2}} \infint e^{-\frac{y_1^2}{2\sigma_2^2}}\ e^{-i y_1 s_1}dy_1 %------------------------------------------------------------------------------- \end{eqnarray}

積分変数$y_1$を含まない係数は,積分の外に出して単純にした. そうすると積分しないといけない,ヤバそうな部分は結構シンプルになる. 式(\ref{eq:8})と(\ref{eq:9})の2つあるけど,よく見ると$\sigma_1$と$\sigma_2$が違うだけで全く同じ形.

積分の部分を計算する

その肝心の積分変数は,「ガウス関数のフーリエ変換はガウス関数になる」という公式が使えて実際には積分計算をする必要がない.

\begin{eqnarray} \mathcal{F}\left[\frac{1}{\sqrt{2\pi}\sigma}e^{-t^2/2\sigma^2}\right](\omega) =\frac{1}{\sqrt{2\pi}\sigma} \infint e^{-t^2/2\sigma^2} e^{-i \omega t}dt =e^{-\sigma^2\omega^2/2} \end{eqnarray}

この公式の変数を置換するだけで,式(\ref{eq:8})と(\ref{eq:9})の積分の部分を計算できる.

\begin{eqnarray} \infint e^{-\frac{y_1^2}{2\sigma_1^2}}\ e^{-i y_1 s_1}dy_1 =\sqrt{2\pi}\sigma_1\ e^{-\frac{\sigma_1^2 s_1^2}{2}}\\ %------------------------------------------------------------------------------- \infint e^{-\frac{y_1^2}{2\sigma_2^2}}\ e^{-i y_1 s_1}dy_1 =\sqrt{2\pi}\sigma_2\ e^{-\frac{\sigma_2^2 s_1^2}{2}} \end{eqnarray}

一重目のフーリエ変換を計算

これで, 式(\ref{eq:6})に上の結果を代入すると, 式(\ref{eq:4})の内側のフーリエ変換はちゃんと計算できる.

\begin{eqnarray} \infint f(y_1, y_2)\ e^{-i y_1 s_1}dy_1 &=& \frac{1}{2\pi\sigma_1^2} e^{-\frac{y_2^2}{2\sigma_1^2}} \ \sqrt{2\pi}\sigma_1\ e^{-\frac{\sigma_1^2 s_1^2}{2}}\\ &-& \frac{1}{2\pi\sigma_1\sigma_2} e^{-\frac{y_2^2}{2\sigma_1^2}} \ \sqrt{2\pi}\sigma_2\ e^{-\frac{\sigma_2^2 s_1^2}{2}}\\ %------------------------------------------------------------------------------- &=& \frac{1}{\sqrt{2\pi}\sigma_1} e^{-\frac{y_2^2}{2\sigma_1^2}}\ e^{-\frac{\sigma_1^2 s_1^2}{2}}\\ &-& \frac{1}{\sqrt{2\pi}\sigma_1} e^{-\frac{y_2^2}{2\sigma_1^2}}\ e^{-\frac{\sigma_2^2 s_1^2}{2}}\\ %------------------------------------------------------------------------------- &=& \label{eq:15} \frac{1}{\sqrt{2\pi}\sigma_1} e^{-\frac{y_2^2}{2\sigma_1^2}} \left( e^{-\frac{\sigma_1^2 s_1^2}{2}}-e^{-\frac{\sigma_2^2 s_1^2}{2}} \right) \end{eqnarray}

二重目のフーリエ変換を計算

ようやく計算し終えた一重目のフーリ変換(式(eq:15))を元の二次元フーリ変換の式(\ref{eq:4})に代入する.

\begin{eqnarray} F(s_1, s_2)&=& \infint \frac{1}{\sqrt{2\pi}\sigma_1}e^{-\frac{y_2^2}{2\sigma_1^2}} \left( e^{-\frac{\sigma_1^2 s_1^2}{2}}-e^{-\frac{\sigma_2^2 s_1^2}{2}} \right) e^{-i y_2 s_2}dy_2\\ %------------------------------------------------------------------------------- &=& \frac{1}{\sqrt{2\pi}\sigma_1} \left( e^{-\frac{\sigma_1^2 s_1^2}{2}}-e^{-\frac{\sigma_2^2 s_1^2}{2}} \right) \infint e^{-\frac{y_2^2}{2\sigma_1^2}} \ e^{-i y_2 s_2}dy_2 \end{eqnarray}

また,積分の部分がガウス関数のフーリ変換の形になっているので,一重目と同様にして公式を使うと計算できる.

\begin{eqnarray} %------------------------------------------------------------------------------- F(s_1, s_2)&=& \frac{1}{\sqrt{2\pi}\sigma_1} \left( e^{-\frac{\sigma_1^2 s_1^2}{2}}-e^{-\frac{\sigma_2^2 s_1^2}{2}} \right) \sqrt{2\pi}\sigma_1 e^{-\frac{\sigma_1^2 s_2^2}{2}}\\ %------------------------------------------------------------------------------- &=& \label{eq:20} e^{-\frac{\sigma_1^2 s_2^2}{2}} \left( e^{-\frac{\sigma_1^2 s_1^2}{2}}-e^{-\frac{\sigma_2^2 s_1^2}{2}} \right) %------------------------------------------------------------------------------- \end{eqnarray}

論文中のフーリエ変換の式と比較

論文[2]の式は次のようになっていて,上の式と微妙に違うので,正しいの確認してみる.

\begin{eqnarray} \label{eq:21} F(s_1, s_2)&=& e^{-2\pi^2\sigma_1^2 t_2^2} \left( e^{-2\pi^2\sigma_1^2 t_1^2}-e^{-2\pi^2\sigma_2^2 t_1^2} \right) \end{eqnarray}

まず$(t_1, t_2)$というのが出てくるけど,実際の論文では回転を考慮しているので,$(s_1, s_2)$を回転させた$(t_1,t_2)$というのが出てくる. 予め回転させた画像を与えると,ここは関係ないので,これは問題ない. あと,$2\pi^2$というなぞの係数が登場するけど,これはフーリエ変換の定義で周波数を$\omega=2\pi k$のように, 単位画素あたり位相$\omega$ではなてく,波数$k$で定義すると現れる係数.要するに単位の違い. 実際に式(\ref{eq:20})に$s_1=2\pi k_1$と$s_2 =2\pi k_2$を代入すると,式(\ref{eq:21})と同じ形になる.

参考文献

  1. Barbara Blakeslee and Mark E. McCourt: "A multiscale spatial filtering account of the White effect, simultaneous brightness contrast and grating induction," Vision Research, vol. 39, issue 26, pp. 4361-4377, 1999.
    doi: 10.1016/S0042-6989(99)00119-4
  2. Barbara Blakeslee, Davis Cope, and Mark E. McCourt1: "The Oriented Difference of Gaussians (ODOG) Model of Brightness Perception: Overview and Executable," Behavior Research Methods, pp. 306–312, vol. 48, no. 1, 2016.
    doi: 10.3758/s13428-015-0573-4
  3. "Difference of Gaussians" in Wikipedia
    https://en.wikipedia.org/wiki/Difference_of_Gaussians