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）に変えると，一連の処理は次のようになる．
- 入力画像$f$に対してFFTすることにより$P=\mathcal{F}[p]$を計算する．
- $F_k = \sum_{l=0}^6 w_l\cdot \mathcal{F}[f_{k,l}]$を計算して配列に蓄積する．
- $\mathcal{F}[Q_k]$を計算する．
- $\mathcal{F}[Q_k]$の全画素の自乗和により$E_k$を計算する．
- $\mathcal{F}[R]$を計算する．
- 逆高速フーリエ変換により$R$を計算する．
参考文献
-
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 - "Difference of Gaussians" in Wikipedia
https://en.wikipedia.org/wiki/Difference_of_Gaussians - 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
