Compute Shader でFFTと畳み込み演算でブラー

はじめに

前回,ComputeShaderを利用してSummed Area Tableを作りましたが, 今回はFFT(高速フーリエ変換)をしてみたいと思います。 フーリエ変換ができれば,大きなカーネルの畳み込み計算を現実的な時間で出来て, Unreal Engineで言うところの,Convolution Bloomのようなエフェクトを実装できます。

https://docs.unrealengine.com/latest/INT/Engine/Rendering/PostProcessEffects/Bloom/

高速フーリエ変換

フーリエ変換の表式はいろいろありますが,今回は連続の空間だと以下のもので表されるとします。

{ \displaystyle
\begin{eqnarray}
F(k) = \int f(x) e^{-ikx} dx
\end{eqnarray}
}

離散空間でも,いろいろ表式はありますが,{N}個のデータがあるとして

{ \displaystyle
\begin{eqnarray}
F(t) = \sum_{x=0}^{N-1} f(x) e^{-i\frac{2\pi x}{N}t}
\end{eqnarray}
}

{N\times N}サイズの画像にフーリエ変換をかける場合は

{ \displaystyle
\begin{eqnarray}
F(s,t) = \sum_{x=0}^{N-1}\sum_{y=0}^{N-1} f(x,y) e^{-i\frac{2\pi x}{N}t} e^{-i\frac{2\pi y}{N}s}
\end{eqnarray}
}

これはX方向に1次元のフーリエ変換を行ってから,Y方向にフーリエ変換を行えばよいことになるのですが,
何も考えずに計算をすると,1ピクセルあたり{O(N^2)}の計算量となります。

これを,フーリエ変換の係数の重みの周期性に着目して加算の仕方を工夫したものがFFT(高速フーリエ変換)です。 1ピクセルあたりの計算量が{O(N\log_2{(N)})}まで減ります。

https://ja.wikipedia.org/wiki/%E9%AB%98%E9%80%9F%E3%83%95%E3%83%BC%E3%83%AA%E3%82%A8%E5%A4%89%E6%8F%9B

FFTの実装

ComputeShaderでFFTを行う実装は,以下のサイトを参考にして実装しました。

https://software.intel.com/en-us/articles/fast-fourier-transform-for-image-processing-in-directx-11

コードはこんな感じです。

前回のSATと同様に,共有メモリをうまく使って高速化をしています。 インテルの実装だと,バタフライ演算でダブルバッファを利用しているので,共有メモリを2倍消費しますが,無駄なバリアが発生せずに高速です。
(何故バタフライ演算のインデックスの計算でreversebits命令を使うとうまくいくのか,あんまりよくわかりません。)

FFTの適用結果

512x512サイズの画像に対して,FFTを適用しました。 FFTを適用すると,実数成分の他に虚数成分も出てくるので,パワースペクトルとして{|F(s,t)|^2}を出力しています。

f:id:hikita12312:20171111163807p:plain:w300
元画像
f:id:hikita12312:20171111164042p:plain:w300
FFT適用後パワースペクトル

512x512サイズのFFTGPU実行時間は,GTX970の環境で0.32ms程度でした。

畳み込み演算

次に畳み込み演算について考えてみます。 畳み込み演算の表式は

{ \displaystyle
\begin{eqnarray}
(f \ast g)(x) = \int f(x) g(x-t) dt
\end{eqnarray}
}

ブルーム表現などで行っていたガウスブラーなどの処理は, この畳み込み演算を離散的に行っていたという事になります。

ところで,畳み込み演算に対してフーリエ変換を行うと,以下のように畳み込む関数をフーリエ変換した関数の積になります。

{ \displaystyle
\begin{eqnarray}
\int (f \ast g)(x) e^{-ikx} dx &=& \int dx \int dt f(x) g(x-t) e^{-ikx} \\\
&=& \int f(x) e^{-ikx} dx \int g(x) e^{-ikx} dx \\\
&=& F(k)G(k)
\end{eqnarray}
}

つまり,入力テクスチャのフーリエ変換結果と,フィルタ画像のフーリエ変換結果を用意して,ピクセル同士の積を行った結果を 逆フーリエ変換することで,畳み込み演算の結果を取得できます。
まともにカーネルを作って畳み込み計算をおこなうと,カーネルサイズに比例して計算速度が下がりますが, フーリエ変換を行えばどのようなサイズのカーネルでも一定の時間で計算が可能で,縮小バッファも用いないのでキレイです。

畳み込みブラーの適用結果

まずはガウスブラーをかけてみます。 計算するとわかりますが,ガウス関数フーリエ変換してもガウス関数です。

{ \displaystyle
\begin{eqnarray}
\frac{1}{\sqrt{2\pi}\sigma} e^{-\frac{x^2}{2\sigma^2}} \rightarrow e^{-\frac{\sigma^2x^2}{2}}
\end{eqnarray}
}

標準偏差{\sigma}ガウスブラーをかけるときは,入力テクスチャをフーリエ変換した結果の中央からの距離{r}の位置に, {e^{-\frac{\sigma^2r^2}{2}}}の重みを掛けてから逆フーリエ変換することで,どのようなサイズのガウスブラーもかけることができます。

f:id:hikita12312:20171111180211p:plain:w300
{\sigma=32}
f:id:hikita12312:20171111180234p:plain:w300
{\sigma=128}

512x512サイズのテクスチャに対して,FFTと畳み込みブラーを実行した結果です。 見やすくなるように,トーンマップとガンマ補正を適用しています。 512x512サイズの畳み込み演算のGPU実行時間は,GTX970の環境で0.13ms程度でした。

大雑把に計算して,FFTに0.3ms,畳み込みに0.13ms,逆FFTで0.3msで0.73ms程度のGPU負荷で,どんなサイズのブラーも掛けられるということになります。 実際は入力のFFT結果はブルーム処理やDoF処理で使いまわすので,実用上では1枚の縮小バッファブラーあたり0.5ms程度となるのでしょうか。 決して速くはありあせんが,綺麗な結果を得られることを考えると,良い方法な気がします。

ついでに,ディスクブラーFFTとの畳込みでやってみました。 円形開口のマスクのフーリエ変換はベッセル関数とかになった記憶があるんですが,面倒なのでこれはそのまま円形マスクテクスチャをFFTした結果を利用します。

f:id:hikita12312:20171111175939p:plain:w300
{r=16}
f:id:hikita12312:20171111175958p:plain:w300
{r=32}

大きなディスクブラーもキレイに出来ています。普通にカーネルを作ってディスクブラーを行おうとすると, 13x13のディスクブラーでも49サンプルぐらいしないとダメだったのでそれを考えるとFFTの方が良いのかもしれません。

また,ディスクブラーやガウスブラーみたいな定まった形ではなく,好きなカーネルを使えるのもメリットです。

f:id:hikita12312:20171111174923p:plain:w128

例えば,星型のカーネルを用意して,このカーネルフーリエ変換したもので畳み込めば...

f:id:hikita12312:20171111175642p:plain:w600

ギザギザとした感じのフィルタがかけられます。 どんなに複雑な形状でも,計算速度は変わりません。

なんだか,学校の課題みたいなネタでしたが,畳み込みブラーをうまく使ったら高品質なブルームやDoFができそうです。