はじめに
前回,ComputeShaderを利用してSummed Area Tableを作りましたが,
今回はFFT(高速フーリエ変換)をしてみたいと思います。
フーリエ変換ができれば,大きなカーネルの畳み込み計算を現実的な時間で出来て,
Unreal Engineで言うところの,Convolution Bloomのようなエフェクトを実装できます。
https://docs.unrealengine.com/latest/INT/Engine/Rendering/PostProcessEffects/Bloom/
フーリエ変換の表式はいろいろありますが,今回は連続の空間だと以下のもので表されるとします。
離散空間でも,いろいろ表式はありますが,個のデータがあるとして
サイズの画像にフーリエ変換をかける場合は
これはX方向に1次元のフーリエ変換を行ってから,Y方向にフーリエ変換を行えばよいことになるのですが,
何も考えずに計算をすると,1ピクセルあたりの計算量となります。
これを,フーリエ変換の係数の重みの周期性に着目して加算の仕方を工夫したものがFFT(高速フーリエ変換)です。
1ピクセルあたりの計算量がまで減ります。
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
ComputeShaderでFFTを行う実装は,以下のサイトを参考にして実装しました。
https://software.intel.com/en-us/articles/fast-fourier-transform-for-image-processing-in-directx-11
コードはこんな感じです。
#if defined(PASS_1ST) && !defined(INVERSE)
#define INPUT_REALPART
#endif
#if defined(PASS_2ND) && defined(INVERSE)
#define OUTPUT_REALPART
#endif
#if LENGTH<=512
#define DOUBLE_BUFFER
#endif
Texture2D<float3> g_SrcTextureR : register(t0);
#if !defined(INPUT_REALPART)
Texture2D<float3> g_SrcTextureI : register(t1);
#endif
RWTexture2D<float3> g_DstTextureR : register(u0);
#if !defined(OUTPUT_REALPART)
RWTexture2D<float3> g_DstTextureI : register(u1);
#endif
static const uint REAL_PART = 0;
static const uint IMAGINARY_PART = 1;
#ifdef DOUBLE_BUFFER
groupshared float3 s_aaTempBuffer[2][2][LENGTH];
#define TempBuffer(tempId, x, ri) (s_aaTempBuffer[(tempId)][(ri)][(x)])
#else
groupshared float3 s_aaTempBuffer[2][LENGTH];
#define TempBuffer(tempId, x, ri) (s_aaTempBuffer[(ri)][(x)])
#endif
void GetButterflyValues(uint passIndex, uint x, out uint2 indices, out float2 weights)
{
int sectionWidth = 2 << passIndex;
int halfSectionWidth = sectionWidth / 2;
int sectionStartOffset = x & ~(sectionWidth - 1);
int halfSectionOffset = x & (halfSectionWidth - 1);
int sectionOffset = x & (sectionWidth - 1);
sincos(2.0f*PI*sectionOffset / (float)sectionWidth, weights.y, weights.x);
weights.y = -weights.y;
indices.x = sectionStartOffset + halfSectionOffset;
indices.y = sectionStartOffset + halfSectionOffset + halfSectionWidth;
if (passIndex == 0)
{
indices = reversebits(indices) >> (32 - BUTTERFLY_COUNT) & (LENGTH - 1);
}
}
void ButterflyPass(int passIndex, uint x, uint tempId, out float3 resultR, out float3 resultI)
{
uint2 indices;
float2 weights;
GetButterflyValues(passIndex, x, indices, weights);
const float3 inputR1 = TempBuffer(tempId, indices.x, REAL_PART);
const float3 inputI1 = TempBuffer(tempId, indices.x, IMAGINARY_PART);
const float3 inputR2 = TempBuffer(tempId, indices.y, REAL_PART);
const float3 inputI2 = TempBuffer(tempId, indices.y, IMAGINARY_PART);
#ifndef DOUBLE_BUFFER
GroupMemoryBarrierWithGroupSync();
#endif
#ifdef INVERSE
resultR = (inputR1 + weights.x * inputR2 + weights.y * inputI2) * 0.5f;
resultI = (inputI1 - weights.y * inputR2 + weights.x * inputI2) * 0.5f;
#else
resultR = inputR1 + weights.x * inputR2 - weights.y * inputI2;
resultI = inputI1 + weights.y * inputR2 + weights.x * inputI2;
#endif
}
#ifdef OUTPUT_REALPART
void InverseButterflyPassRealPart(int passIndex, uint x, uint tempId, out float3 resultR)
{
uint2 indices;
float2 weights;
GetButterflyValues(passIndex, x, indices, weights);
const float3 inputR1 = TempBuffer(tempId, indices.x, REAL_PART);
const float3 inputR2 = TempBuffer(tempId, indices.y, REAL_PART);
const float3 inputI2 = TempBuffer(tempId, indices.y, IMAGINARY_PART);
resultR = (inputR1 + weights.x * inputR2 + weights.y * inputI2) * 0.5f;
}
#endif
[numthreads(LENGTH, 1, 1)]
void CS_Transform(uint3 position : SV_DispatchThreadID)
{
const uint bufferId = position.x;
#ifdef PASS_1ST
uint2 texturePos = position.xy;
#endif
#ifdef PASS_2ND
uint2 texturePos = position.yx;
#endif
texturePos = (texturePos + LENGTH / 2) % LENGTH;
TempBuffer(0, bufferId, REAL_PART) = g_SrcTextureR[texturePos];
#if defined(INPUT_REALPART)
TempBuffer(0, bufferId, IMAGINARY_PART) = 0.0f;
#else
TempBuffer(0, bufferId, IMAGINARY_PART) = g_SrcTextureI[texturePos];
#endif
for (uint i = 0; i < BUTTERFLY_COUNT - 1; i++)
{
GroupMemoryBarrierWithGroupSync();
const uint srcTempId = i % 2;
const uint dstTempId = (i+1) % 2;
ButterflyPass(i, bufferId, srcTempId, TempBuffer(dstTempId, bufferId, REAL_PART), TempBuffer(dstTempId, bufferId, IMAGINARY_PART));
}
{
GroupMemoryBarrierWithGroupSync();
const uint srcTempId = (BUTTERFLY_COUNT - 1)% 2;
#ifdef OUTPUT_REALPART
InverseButterflyPassRealPart(BUTTERFLY_COUNT - 1, bufferId, srcTempId, g_DstTextureR[texturePos]);
#else
ButterflyPass(BUTTERFLY_COUNT - 1, bufferId, srcTempId, g_DstTextureR[texturePos], g_DstTextureI[texturePos]);
#endif
}
}
前回のSATと同様に,共有メモリをうまく使って高速化をしています。
インテルの実装だと,バタフライ演算でダブルバッファを利用しているので,共有メモリを2倍消費しますが,無駄なバリアが発生せずに高速です。
(何故バタフライ演算のインデックスの計算でreversebits命令を使うとうまくいくのか,あんまりよくわかりません。)
512x512サイズの画像に対して,FFTを適用しました。
FFTを適用すると,実数成分の他に虚数成分も出てくるので,パワースペクトルとしてを出力しています。
512x512サイズのFFTのGPU実行時間は,GTX970の環境で0.32ms程度でした。
畳み込み演算
次に畳み込み演算について考えてみます。
畳み込み演算の表式は
ブルーム表現などで行っていたガウスブラーなどの処理は,
この畳み込み演算を離散的に行っていたという事になります。
ところで,畳み込み演算に対してフーリエ変換を行うと,以下のように畳み込む関数をフーリエ変換した関数の積になります。
つまり,入力テクスチャのフーリエ変換結果と,フィルタ画像のフーリエ変換結果を用意して,ピクセル同士の積を行った結果を
逆フーリエ変換することで,畳み込み演算の結果を取得できます。
まともにカーネルを作って畳み込み計算をおこなうと,カーネルサイズに比例して計算速度が下がりますが,
フーリエ変換を行えばどのようなサイズのカーネルでも一定の時間で計算が可能で,縮小バッファも用いないのでキレイです。
畳み込みブラーの適用結果
まずはガウスブラーをかけてみます。
計算するとわかりますが,ガウス関数はフーリエ変換してもガウス関数です。
標準偏差がのガウスブラーをかけるときは,入力テクスチャをフーリエ変換した結果の中央からの距離の位置に,
の重みを掛けてから逆フーリエ変換することで,どのようなサイズのガウスブラーもかけることができます。
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した結果を利用します。
大きなディスクブラーもキレイに出来ています。普通にカーネルを作ってディスクブラーを行おうとすると,
13x13のディスクブラーでも49サンプルぐらいしないとダメだったのでそれを考えるとFFTの方が良いのかもしれません。
また,ディスクブラーやガウスブラーみたいな定まった形ではなく,好きなカーネルを使えるのもメリットです。
例えば,星型のカーネルを用意して,このカーネルをフーリエ変換したもので畳み込めば...
ギザギザとした感じのフィルタがかけられます。
どんなに複雑な形状でも,計算速度は変わりません。
なんだか,学校の課題みたいなネタでしたが,畳み込みブラーをうまく使ったら高品質なブルームやDoFができそうです。