モーションブラーの実装

モーションブラー

ゲームのフレームとフレームの間の動きに応じて,ブラーを与える効果がモーションブラーです。 レンダリングされた一枚のシーンがカメラによって撮影されたものと考えて,露出時間からブラーを計算します。

https://docs.unity3d.com/ja/current/Manual/PostProcessing-MotionBlur.html

モーションブラーを適用することで,レースゲームのスピード感が出たり, フレームレートが低い場合でも,滑らかに動いているような表現ができます。

今回は,モーションブラーの方法として,「A Reconstruction Filter for Plausible Motion Blur」という論文を内容を実装してみます。

http://casual-effects.com/research/McGuire2012Blur/index.html

「A Reconstruction Filter for Plausible Motion Blur」の方法はUnityでも実装されている方法のようです。

速度マップ

モーションブラーを作成するために,オブジェクトがどのような速度を持っているかの情報を持った速度マップを作成します。
一般に３Dモデルを画面に表示するときには,Projection,View,Modelといった,3つの行列を利用して,ワールド空間からスクリーン空間へと変換することで実現します。 このとき1フレーム前の行列を保存しておいて,現在のフレームとの差分を考えることで,フレーム間のオブジェクトの速度が計算されます。

float4x4 g_CurrentModelMatrix, g_CurrentViewMatrix, g_CurrentProjectionMatrix;
float4x4 g_PrevModelMatrix, g_PrevViewMatrix, g_PrevProjectionMatrix;

float4 currentPosition = position * CurrentModelMatrix * CurrentViewMatrix * CurrentProjectionMatrix;
float4 prevPosition = position * PrevModelMatrix * PrevViewMatrix * PrevProjectionMatrix;

float2 velocity = currentPosition.xy/currentPosition.w - prevPosition.xy/prevPosition.w;
velocity *= 0.5f; // 画面サイズを1.0としたUV座標系での速度に変換

画像は右に移動をするオブジェクトに対して,速度マップを作成した結果です。 見やすくするために,実際の速度の20倍の値を表示しています。

A Reconstruction Filter for Plausible Motion Blur

先程の速度マップの速度情報に従ってボカす方向とボカす量を決めればモーションブラーが出来上がります。 しかし,速度マップは,そのフレームでオブジェクトがあった位置にしか,速度情報は書き込まれません。

上の図のような状態では,濃い青色のオブジェクトの部分には「右に移動している」という情報が格納されていますが, 本来通過してきたはずのオブジェクトの左側の境界の外の付近の速度マップの値は0です。 単純に現在の速度マップの値を信用してブラーをかけると,境界外部に対してブラーをかけることが出来ません。

A Reconstruction Filter for Plausible Motion Blurの方法では, 注目するピクセルの周辺を含めた支配的な速度を仮定してブラーを行って,境界外部へのブラーを解決しています。 紹介されている方法では3つのパスを利用してモーションブラーを適用します。 まず,サイズの速度マップがあったとして,サイズのTileMaxと呼ばれるバッファを作成します。

TileMaxは,元の速度マップをのタイルに分割して,そのタイル内部で最大の速度を格納します。 図の赤い矢印が中心のタイルが保持する最大の速度です。
さらに,TileMaxから周辺のタイルを含めて最大の速度を持ったTileを自身に格納する,サイズのNeighborMaxと呼ばれるバッファも作成します。

図の赤い矢印が,周辺を含めた最大の速度として,NeighborMaxになります。
この,NeighborMaxがブラーを行うときに,支配的な周辺の速度であるとして,利用されます。

最後に,シーンのカラーテクスチャ,深度テクスチャ,元の速度マップ,そして今作成されたNeighborMaxを利用して,最終的なブラー結果を作成します。 シェーダーは論文のコードをそのまま実装しただけですが,以下のようなものになります。

クリックで展開

Texture2D g_ColorTexture : register(t0);
Texture2D g_VelocityTexture : register(t1);
Texture2D g_DepthTexture : register(t2);
Texture2D g_NeighborMax : register(t3);
SamplerState g_LinearSampler : register(s0);

#define    PIXEL_SIZE g_aTemp[0].xy
#define    SAMPLE_NUM g_aTemp[0].z
#define    BLUR_SCALE g_aTemp[0].w
#define    K g_aTemp[1].x
#define    VELOCITY_SCALE g_aTemp[1].y
#define    DEPTH_FACTOR_1      g_aTemp[1].z
#define    DEPTH_FACTOR_2      g_aTemp[1].w

static const float EPSILON = 0.0001f;
static const float SOFT_Z_EXTENT = 0.1f;

float GetLinearDepth(float2 uv)
{
    const float depth = SampleLod0(g_DepthTexture, g_LinearSampler, uv).r;
    const float distance = DepthToViewPositionZ(depth, DEPTH_FACTOR_1, DEPTH_FACTOR_2);
    return -distance; // 論文だとz-が奥方向
}

float2 GetVelocity(Texture2D velocityMap, float2 uv)
{
    float2 velocity = SampleLod0(g_NeighborMax, g_LinearSampler, uv).rg; // 画面を1.0としたUV座標系Velocity
    velocity *= VELOCITY_SCALE; // スケーリング
    const float v = length(velocity);
    velocity = velocity / (v + EPSILON) * clamp(v, 0.5f*PIXEL_SIZE.y, K*PIXEL_SIZE.y); // 0.5px~Kpxのclamp
    return velocity;
}

float softDepthCompare(float za, float zb)
{
    return clamp(1.0f - (za - zb) / SOFT_Z_EXTENT, 0.0f, 1.0f);
}

float cone(float L, float mag)
{
    return clamp(1.0f - L / mag, 0.0f, 1.0f);
}

float cylinder(float L, float mag)
{
    return 1.0f - smoothstep(0.95f*mag, 1.05f*mag, L);
}

float4 PS_Reconstruct(PPFX_OUT In) : SV_TARGET0
{
    const float2 X = In.UV.xy;
    const float3 color = SampleLod0(g_ColorTexture, g_LinearSampler, X).rgb;
    const float2 neighborMaxVelocity = GetVelocity(g_NeighborMax, X);
    // 1px以下のブラーならば,何もしない
    const float2 harfPixelSize = PIXEL_SIZE * 0.5f;
    if (length(neighborMaxVelocity) < EPSILON + harfPixelSize.y)
    {
        return float4(color.r, color.g, color.b, 1.0f);
    }
    // Sample the current point
    const float2 velocityX = GetVelocity(g_VelocityTexture, X);
    const float velocityXLength = length(velocityX);
    float weight = 1.0f / velocityXLength;
    float3 sum = color.rgb * weight;
    // Take S-1 additional neighbor samples
    const float j = rand(X) - 0.5f; // jitter random
    const int selfIndex = ((int)(SAMPLE_NUM) - 1) / 2;
    for (int i = 0; i < SAMPLE_NUM; ++i)
    {
        if (i == selfIndex) continue;
        const float t = lerp(-1.0f, 1.0f, ((float)(i)+j+1.0f)/(SAMPLE_NUM+1.0f));
        // round to nearest
        const float2 Y = X + neighborMaxVelocity * t + harfPixelSize;
        // Fore- vs. background classification of Y relative to X
        const float zx = GetLinearDepth(X);
        const float zy = GetLinearDepth(Y);
        const float2 velocityY = GetVelocity(g_VelocityTexture, Y);
        const float velocityYLength = length(velocityY);
        const float L = length(X - Y);
        float alphaY =
            // Case 1: Blurry Y in front of any X
            softDepthCompare(zx, zy) * cone(L, velocityYLength)
            // Case 2: Any Y behind blurry X; estimate background
            + softDepthCompare(zy, zx) * cone(L, velocityXLength)
            // Case 3: Simultaneously blurry X and Y
            + cylinder(L, velocityYLength) * cylinder(L, velocityXLength) * 2.0f
            ;
        alphaY *= BLUR_SCALE;
        // Accumulate
        weight += alphaY;
        sum += alphaY * SampleLod0(g_ColorTexture, g_LinearSampler, Y).rgb;
    }
    sum = sum / weight;
    return float4(sum.r, sum.g, sum.b, 1.0f);
}

シェーダーでは,NeighborMaxの速度ベクトル方向に対してサンプリングを行っていき,実際のサンプリングされた地点の速度と深度値からブラー量を決定しています。 この時,注目しているテクセルよりも,サンプルされたテクセルが前にあるか後にあるかでどの重みを採用するか計算を行っているようです。 また,jitterとして,ランダムなオフセットを与えることで,サンプル数が少ない場合でもそれなりに見えるようになっています。

論文ではOpenGL座標系だからか,カメラ奥方向がZ-となっていたので修正しています。 また,ブラー量を無理やりスケールする係数BLUR_SCALEも追加しています。

結果

結果です。BLUR_SCALE=6,K=16,SAMPLE_NUM=16で適用しました。 GPU実行時間は,GTX970の環境で1280x720サイズのテクスチャで.0.25ms程度でした。

Motion Blur OFF	Motion Blur ON

Motion Blur OFF	Motion Blur ON

まだまだキレイにする工夫の余地はありそうですが,とりあえずこんな感じです。

FXAA3の組み込み

FXAA

スクリーンスペースのアンチエイリアスの処理として,NVIDIAのFXAAが有名です。
ソースはNVIDIAのGameWorksのサンプルにあります。

https://developer.nvidia.com/gameworks-directx-samples

現在FXAAのバージョンは3.11のようで,FXAAの3系はマクロによって,スケーラブルにクォリティを決定できます。
細かい実装は把握していませんが,テクスチャの輝度情報からエッジを検出し,エッジ部分をボカすような事を行っているような気がします。

何よりも,1パスで手軽にそれなりのアンチエイリアスが掛けられるのが利点かと思います。 ちなみに,UnityのPost Processing StackもFXAAでアンチエイリアスを掛けているようです。

FXAAの組み込み

まず,NVIDIAのGameWorksを落としてきて,FXAA.hlslのサンプルコードを持ってきます。 FXAA.hlslのコメントにかかれている通りに組み込めば良いだけですが輝度について,少しだけ修正を行いました。

FXAA.hlslの725行目付近に以下のような定義があります。

#if (FXAA_GREEN_AS_LUMA == 0)
    FxaaFloat FxaaLuma(FxaaFloat4 rgba) { return rgba.w; }
#else
    FxaaFloat FxaaLuma(FxaaFloat4 rgba) { return rgba.y; }
#endif   

FXAA_GREEN_AS_LUMAマクロによって,テクスチャのどの成分を輝度値とみなすかを切り替える関数ですが, 輝度値をAlpha成分に格納したり,G成分で代用するよりも,直接輝度の計算を行うべきかと思ったので,内積計算に置き換えます。

FxaaFloat FxaaLuma(FxaaFloat4 rgba) { return dot(rgba.rgb, float3(0.298912f, 0.586611f, 0.114478f)); }

あとは,以下のようにピクセルシェーダーのパスに組み込めば,FXAAを利用する準備ができました。

#define FXAA_PC 1
#define FXAA_HLSL_5 1
//#define FXAA_QUALITY__PRESET 12
#include "FXAA.hlsl"

#define    PIXEL_SIZE g_aTemp[0].xy
#define    SUBPIX g_aTemp[1].x
#define    EDGE_THRESHOLD g_aTemp[1].y
#define    EDGE_THRESHOLD_MIN g_aTemp[1].z

float4 PS_FXAA3(PPFX_OUT In) : SV_TARGET0
{
    float2 uv = In.UV.xy;
    FxaaTex InputFXAATex = { g_LinearSampler, g_SrcTexture };
    return FxaaPixelShader(
        uv,                         // FxaaFloat2 pos,
        FxaaFloat4(0.0f, 0.0f, 0.0f, 0.0f),     // FxaaFloat4 fxaaConsolePosPos,
        InputFXAATex,                   // FxaaTex tex,
        InputFXAATex,                   // FxaaTex fxaaConsole360TexExpBiasNegOne,
        InputFXAATex,                   // FxaaTex fxaaConsole360TexExpBiasNegTwo,
        PIXEL_SIZE,                 // FxaaFloat2 fxaaQualityRcpFrame,
        FxaaFloat4(0.0f, 0.0f, 0.0f, 0.0f),     // FxaaFloat4 fxaaConsoleRcpFrameOpt,
        FxaaFloat4(0.0f, 0.0f, 0.0f, 0.0f),     // FxaaFloat4 fxaaConsoleRcpFrameOpt2,
        FxaaFloat4(0.0f, 0.0f, 0.0f, 0.0f),     // FxaaFloat4 fxaaConsole360RcpFrameOpt2,
        SUBPIX,                     // FxaaFloat fxaaQualitySubpix,
        EDGE_THRESHOLD,             // FxaaFloat fxaaQualityEdgeThreshold,
        EDGE_THRESHOLD_MIN,         // FxaaFloat fxaaQualityEdgeThresholdMin,
        0.0f,                          // FxaaFloat fxaaConsoleEdgeSharpness,
        0.0f,                          // FxaaFloat fxaaConsoleEdgeThreshold,
        0.0f,                          // FxaaFloat fxaaConsoleEdgeThresholdMin,
        FxaaFloat4(0.0f, 0.0f, 0.0f, 0.0f)      // FxaaFloat fxaaConsole360ConstDir,
    );
}

PS3とか,Xbox360のコンソール機向けのパラメータがありますが,基本的に0.0で問題ないと思います。

結果

FXAA_QUALITY__PRESET=12として,FXAAをかけた結果です。

FXAA OFF	FXAA ON

拡大すると分かりやすい気がします。

GPU実行時間は,GTX970の環境で1280x720サイズのテクスチャで.0.15ms程度でした。 FXAA_QUALITY__PRESETを高品質のもの(39とか)にすると,さらに負荷はかかりますが,少しキレイになります。

レンズゴーストの実装

レンズゴースト

レンズゴーストを実装します。
レンズゴーストはブルームと違って,レンズ内部で反射を繰り返した光がハッキリと像となって写り込んだものです。

http://av.jpn.support.panasonic.com/support/dsc/knowhow/knowhow15.html

今回は以下のサイトの実装を参考にしました。

http://john-chapman-graphics.blogspot.jp/2013/02/pseudo-lens-flare.html

実装

実装のアイディア自体は非常に簡単です。

光源となるテクスチャを用意して,画面中心に対して点対称となるようにテクスチャをフェッチすれば,レンズゴーストになります。 このときに,フェッチする回数を増やしたり,カラースケールを変更することでレンズゴーストっぽくなります。 HLSLのコードは以下のような感じになりました。

Texture2D g_SrcTexture : register(t0);
SamplerState g_LinearSampler : register(s0);

#define DISPERSION_SCALE g_aTemp[0].x
#define DISPERSION_OFFSET g_aTemp[0].y
#define GHOST_MASK_RADIUS g_aTemp[0].z
#define GHOST_MASK_FADE g_aTemp[0].w
#define GHOST_COLOR(index) g_aTemp[(index)%4+1].xyz
static const float2 CENTER = float2(0.5f, 0.5f);

float4 PS_Create(PPFX_OUT In) : SV_TARGET0
{
    float4 outColor = float4(0.0f, 0.0f, 0.0f, 1.0f);
    const float2 uv = In.UV.xy;
    const float2 ghostBaseUV = float2(1.0f, 1.0f) - uv;
    const float2 toCenterVec = (CENTER - ghostBaseUV) * DISPERSION_SCALE;
    for (uint index = 0; index < ITERATION_NUM; ++index)
    {
        const float2 ghostUV = ghostBaseUV + toCenterVec * (DISPERSION_OFFSET + index);
        float3 ghostColor = SampleLod0(g_SrcTexture, g_LinearSampler, ghostUV).rgb;

        const float r = length(CENTER - ghostUV);
        const float fadeoutR = (r - GHOST_MASK_RADIUS);
        const float weight =
            (r / length(CENTER)) // 光源が重なったときに極端に明るくさせない重み
            * (1.0f - smoothstep(0.0f, GHOST_MASK_FADE, r - GHOST_MASK_RADIUS)) // 外周をフェードアウトさせる重み
            ;
        outColor.rgb += ghostColor * GHOST_COLOR(index) * weight;
    }

    return outColor;
}

フェッチするマスクのフェードアウト項

     const float weight =
            (r / length(CENTER)) // 光源が重なったときに極端に明るくさせない重み
            * (1.0f - smoothstep(0.0f, GHOST_MASK_FADE, r - GHOST_MASK_RADIUS)) // 外周をフェードアウトさせる重み
            ;

画面の端のほうもフェッチすると,画面の境界で見切れてしまった像が映ってしまうので,端に行くほど滑らかにフェードアウトさせる重みを与えます。 また,画面中心に光源が存在する場合は,中央がとても明るくなってしまうので,適当に減衰させています。

パスの繰り返し
このシェーダーのITERATION_NUMを増やせばゴーストの像はたくさん作られますが,数十回もループさせると負荷が心配です。 そこで,このシェーダーの適用結果を入力にして,繰り返して適用することで,ピンポンブラー的にフェッチ数を増やすことが出来ます。 ITERATION_NUM=8のパスを2回繰り返せば64個のゴーストが生成されます。 また,繰り返すことでゴーストの色合いを作るマスクも繰り返し適用されるので,複雑な色合いも表現できます。

適用結果

高速化の為に,ブルーム効果のフィルタ1レベル目のブラー結果を入力として,ゴースト用のテクスチャは1/4サイズの縮小バッファを利用しています。 像がハッキリと出るタイプのゴーストではないですが,まあまあキレイに適用できている気がします。

シェーダー管理

シェーダー管理

今までいろいろとシェーダーを書いてきましたが,現在行っている管理の仕方を紹介します。

https://www.slideshare.net/siliconstudio/cedec2005-kawase
http://research.tri-ace.com/

Shader Package

前提としてDirectX11環境でHLSLでShader Model 5.0のターゲットでシェーダーを作成しています。 ゲームでシェーダーを利用するときには,コンパイル済みのバイナリを読み込むことになります。 このとき,特定のシーンで利用するシェーダーの名称をKeyとしてシェーダーバイナリを引っ張る辞書的な構造が必要となるのですが, これをShader Pacakgeと呼んで事前に作成しました。

ShaderPackageの構成要素は次の通りです。

Effect

ShaderPackageの大本となる一つのファイルを指します。1ShaderPackageにつき1Effectです。

Technique

1つのシェーダーファイルに対応する構造です。1つのEffectに複数のTechniqueが含まれます。

Pass

シェーダーのマクロ分岐や引数の違いによるバリエーションです。1つのTechniqueに複数のPassが含まれます。

かつてDirextXにあった,fx拡張子のエフェクトシステムに似た構造を利用します。 Passが最小単位で,利用するVertexShader,PixelShaderのエントリポイント名や,マクロ定義などを含みます。

例えば,上の画像のようなEffectファイルを用意したとき,Tonemap.TypeAのPixelShaderは次のように事前にシェーダーコンパイルされることになります。

//// ShaderTonemap.fx
// Vertex Shader Entry
PPFX_OUT VS_Main(PPFX_IN In)
{
}
// Pixel Shader Entry
float4 PS_Main(PPFX_Out In)
{
#if TONEMAP==A
  return FuncA(In);
#endif
#if TONEMAP==B
  return FuncB(In);
#endif
}

fxc.exe ShaderTonemap.fx /Gfa /T ps_5_0 /E PS_Main /D TONEMAP=A /D GAMMA /Fo Tonemap_TypeA_PS.o

ビルドされたシェーダーバイナリを一つのファイルにまとめて, Technique名とPass名をKeyとして,シェーダーバイナリをValueとするような構造としたバイナリデータをShader Pacakgeとして出力します

あとは,実際にゲームでシェーダーを利用するときに,Shader Pacakgeを利用して, 必要となるTechniqueとPassの辞書から必要なシェーダーバイナリを取り出して,DirextX11のAPIに渡してあげればOKです。

シェーダーの組合せ爆発の問題

高速に動作するシェーダーを記述する場合,極力if文による条件分岐やfor文によるイテレーションを避け,シェーダーアセンブリが展開されやすいコードを書くことになります。 そのため,シェーダーのバリエーションを作る場合には,条件分岐ごとにシェーダーを用意する必要があります。 組合せ爆発とは,そのバリエーションが簡単に膨れ上がってしまうことを指します

例えば,トーンマップを例に取って考えてみると,

• トーンマップ関数の種類 5種類 (Linear,Reinhard,exp,log,Filmic)
• カラーコレクション 3種 (OFF, ColorMatrix, LUT)
• Noise ON/OFF 2種
• グレアの合成 ON/OFF 2種
• ガンマ変換 ON/OFF 2種
  (HDRディスプレイ出力とか,色空間変換をすると,もっと増える)

これだけで,5x3x2x2x2=120種類ものシェーダーをビルドして,バイナリに含めなければなりません。 もはやトーンマップのケースで言えば,これは半分しょうがない問題なのかもしれません。

しかし,実際のシェーダーバリエーションとしては,要らない組み合わせのものが含まれたりするかもしれません。 例えば,ゲーム内でトーンマップ関数のうち、Filmic以外はグレア処理は使わないとか....

シェーダー定義の設定ファイルとしてのLua

最初はShaderPacakgeのTechnique,Passの定義ファイルをJSON形式で記述していたのですが,手作業でバリエーションを書くのは辛い作業でした。
そこで,スクリプト言語であるLuaをC++に組み込んで,シェーダーの定義をプログラマブルに行えるようにしました。

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

local Util = require("script.Utility")
root = {}
-- ShaderPackage Name
root["name"] = "ShaderPackage_PFX"
--  Technique
root["technique"] = {
    -- DoFのソース作成パス
    ["DofSource"] = {
        ["src"]  = "src/ppfx/PfxDofSource.fx",
        ["pass"] = {
            ["ColorRGB_CoCA"] = {["vs"]="VS_ScreenRect", ["ps"]="PS_Source"},
        }
    },
    -- トーンマップ
    ["Tonemap"] = {
        ["src"]  = "src/ppfx/PfxTonemap.fx",
        ["pass"] = (function ()
            local pass = {}
            pass = Util:PassCombination(pass, {
                ["$name"]   = "{$1}_{$2}{$3}{$4}",
                ["vs"]      = "VS_ScreenRect",
                ["ps"]      = "PS_Tonemap",
                ["$map"] = {
                    [1] = {
                        {"Exp","/D TONEMAP_FUNC_TYPE=TMF_EXP"},
                        {"Log","/D TONEMAP_FUNC_TYPE=TMF_LOG"},
                        {"Reinhard","/D TONEMAP_FUNC_TYPE=TMF_REINHARD"},
                        {"ACESFilmic","/D TONEMAP_FUNC_TYPE=TMF_ACES"},
                    },
                    [2] = {
                        {"RGB",""},
                        {"LUM","/D LUMBASE"},
                    },
                    [3] = {
                        {"",""},
                        {"_Noise","/D NOISE"},
                    },
                    [4] = {
                        {"",""},
                        {"_CMAT","/D COLOR_CORRECTION"},
                        {"_CMATLUT","/D COLOR_CORRECTION /D COLOR_LUT"},
                    }
                }
            })
            pass["Linear"] = {["vs"]="VS_ScreenRect", ["ps"]="PS_Tonemap", ["define"]="/D TONEMAP_FUNC_TYPE=TMF_LINEAR"}
            pass["Linear_NOISE"] = {["vs"]="VS_ScreenRect", ["ps"]="PS_Tonemap", ["define"]="/D TONEMAP_FUNC_TYPE=TMF_LINEAR /D NOISE"}
            return pass
        end)()
    },
-- ...
}

Util::Passcombinationは,Pass名用のタグとマクロ定義の組から,組み合わせを作成するヘルパー関数です。 単にJSONやXMLのようにデータだけを記述するよりも,自由度が高く,かつ手作業での定義も可能で良い方法だと思います。

あと,JSON形式では使えないコメント行を作れるのも利点です。

Screen Space Ambient Occlusion : Alchemy AO

Ambient Occlusion

アンビエントオクルージョンをやってみます。

http://frederikaalund.com/a-comparative-study-of-screen-space-ambient-occlusion-methods/
http://ambientocclusion.hatenablog.com/entry/2013/11/07/152755

アンビエントオクルージョンは下の図のように,環境光と呼ばれる全ての方向から一様に降り注ぐ光が 奥まった場所や凹んだ部分では物体に遮蔽されて,暗く見えるような効果のことを呼びます。
(詳しい説明は参考URLのページに任せます)

左側の凹んだ場所にある点は,右側の点とくらべて,周りの物体に光が遮蔽されて,多少暗く見えるはずです。 暗くなる度合いは,点の接平面の法線方向の半球面上で,開けている領域がどの程度存在するかの積分計算で求めることが出来ます。

この計算を3D空間で真面目にやらないで,カメラに射影された2Dの空間上でDepthテクスチャや法線テクスチャをうまく使って高速に計算する方法が Screen Space Ambient Occlusion (以下,SSAO)です。

Alchemy AO

真面目に半球面上の積分計算を行うと,とてもリアルタイムでは出来ないので,積分計算を高速に行う近似的な方法がいくつかあるのですが, 今回はAlchemy AOという手法を実装してみます。

http://casual-effects.com/research/McGuire2011AlchemyAO/index.html

Alchemy AOは注目する画素周辺にスクリーンスペースでランダムな点を打って,積分計算を近似する手法です。 ちなみに,Unity標準のポストエフェクトのアセットのAmbient Occlusionは,ソースコードを見る限りAlchemyAOっぽいです。

https://github.com/Unity-Technologies/PostProcessing

具体的な表式を見てみます。

注目するテクセルの法線と,その周辺でランダムにサンプルされた方向ベクトルとの内積の総和を取ることが肝のようです。 max演算によって,半球の裏側はサンプルされないということになっています。 は正規化時に発散しないためのパラメータでしょうか。はよくわかりません。

をランダムにサンプルするためのサンプリング半径は透視変換によって決まります。 例えばDirextX系で,View空間座標をクリッピング空間座標に移すとして,

これより,をView空間でのサンプリング半径と置き換えれば,カメラから遠い場所ほどサンプリング半径が小さくなるよう決められます。

実際にAlchemy AOを実装してみます。

元画像	元画像法線

この元画像に従って,のパラメータでAO値を計算した結果です。 ちなみに1280x720サイズの元画像に対して,1/2サイズの640x360サイズでAOを作っています。

溝っぽい場所が暗くなっています。 GPU実行時間は,GTX970の環境で0.51ms程度でした。

Bilateral Blur

上の画像のAO結果にはノイズが残っているので,そのまま元の画像に乗算合成してしまうと汚い画面になってしまいます。 単純にキレイにするならば,サンプリング数を増やせば良いのですが16サンプル程度でも0.5msかかっている以上,サンプリング数を増やすのは難しそうです。 そこで,ノイズ除去の為にブラーを適用します。

単にガウスブラーを適用すると,ノイズは除去されるのですが,下の画像のように画面全体がボケてしまします。

三次元の構造から、エッジを保ったままノイズを除去するために,Bilateral Blurを利用します。

Bilateral Blurの基本的な表式は,1次元の場合は以下のようなものです

位置の色に対して,標準偏差を利用して畳み込みを行うことでブラーをかけた色を得ます。 は正規化のための分配関数です。 重みの項が存在しなければ,ガウスブラーと一緒ですが,このを工夫することで,エッジを保ったままブラーを適用します。

今回利用したBilateral Blurの重みも加えた計算式は以下のようなものです。

法線とDepth値を利用することで,急激に法線や深度が変化する部分をサンプルしないようにしています。 の項は法線が似た方向を持っていなければ0になるような値で, の項はDepthが近い値ほど大きな値になる重みです。Depth値は線形化したものではなく,非線形なDepthをそのまま利用しています。

この重みの付け方は以下の資料のBirateral Upsamplingの重みを参考にしました

http://developer.amd.com/wordpress/media/2012/10/ShopfMixedResolutionRendering.pdf

を利用して,実際にBilateral Blurをかけた結果です。

Gauss Blur	Bilateral Blur

エッジを残しつつ,ノイズを除去することができました。
GPU実行時間は,GTX970の環境で0.2ms程度でした。

ブラーをかけた結果を元の画像に乗算合成した結果が以下のものです。

元画像	SSAO後

ウサギとウサギの間のくぼんでいるところがちゃんと暗くなっています。 画面左側の金色のウサギの隣のウサギとの間の暗がりがわかりやすいです。 全部合わせてコミコミで1280x720サイズのテクスチャの1/2バッファで1.0ms以下の負荷でAOがかけられました。

本当はBilateral Blurの重みの参考にした資料で説明されている, Bilateral Upsamplingと呼ばれるアップサンプリングで 1/2の縮小バッファで作成されたSSAO結果を合成するべきなのですが, ただのBilinearフィルタでのアップサンプルでも, そんなにアーティファクトが目立たないような気がしたので未実装です。

Compute Shader での実装

今までの実装はPixel Shaderを利用したものでしたが,試しにCompute Shaderでも実装してみました。

AlchemyAOを実装する際に,サンプリング時に毎回Depthテクスチャのフェッチが必要になるのですが, ComputeShader版では,32x32のテクセルにスレッドを割当て,32x32サイズのブロックの周囲に64x64のブロックを4つ作って,Depth値を先にキャッシュしています。 実質4テクセルフェッチで半径48ピクセル範囲の計算を行えるようにすることで,高速化を図っています。

64サンプリングでPixelShader版とComputeShader版を比べてみると,1280x720サイズの1/2縮小バッファで,

• PS版 : 2.03ms
• CS版 : 0.47ms

という結果でした。Compute Shader版のほうが4倍近く速いです。 しかし,キャッシュを取る関係上,キャッシュのサイズ以上のサンプリング半径の部分は不正確になり,ブロックノイズが生じている気がします。 もう少し調整の余地がありそうです。

もろもろのソースコードは以下。

クリックで展開

Texture2D g_DepthTexture : register(t0);
Texture2D g_NormalTexture : register(t1);
SamplerState g_LinearSampler : register(s0);

#define    SAMPLE_NUM g_aTemp[0].x        ///< @brief サンプリング数
#define    SAMPLE_RADIUS g_aTemp[0].y     ///< @brief サンプリング半径(ワールド空間)
#define    INTENSITY g_aTemp[0].z         ///< @brief 強度
#define    EPSILON g_aTemp[0].w           ///< @brief AlchemyAO パラメータ
#define    BIAS g_aTemp[1].x              ///< @brief AlchemyAO パラメータ
#define    DEPTH_FACTOR_1 g_aTemp[1].y    ///< @brief プロジェクション行列,パラメータ proj23 or proj32
#define    DEPTH_FACTOR_2 g_aTemp[1].z    ///< @brief プロジェクション行列,パラメータ proj22
#define    DEPTH_FACTOR_3 g_aTemp[1].w    ///< @brief プロジェクション行列,パラメータ proj00
#define    DEPTH_FACTOR_4 g_aTemp[2].x    ///< @brief プロジェクション行列,パラメータ proj11
#define    ASPECT_RATIO g_aTemp[2].y      ///< @brief プロジェクション行列,アスペクト比
#define    NEAR_Z g_aTemp[2].z            ///< @brief プロジェクション行列,Near
#define    FAR_Z g_aTemp[2].w             ///< @brief プロジェクション行列,Far
#define    TEXTURE_SIZE g_aTemp[3].xy     ///< @brief (CSのみ) テクスチャサイズ
#define    TEXTURE_SIZE_X g_aTemp[3].x        ///< @brief (CSのみ) テクスチャサイズ
#define    TEXTURE_SIZE_Y g_aTemp[3].y        ///< @brief (CSのみ) テクスチャサイズ

#ifdef USE_COMPUTE_SHADER
RWTexture2D<float> g_DstTexture : register(u0);
static const uint BLOCK_SIZE = 32;
/// @brief 32x32サイズのブロックを2x2並べたキャッシュ
groupshared float s_aDepthCache[BLOCK_SIZE * BLOCK_SIZE * 4];

/// @brief 共有メモリにView空間座標を読み込み
void LoadDepthCache(uint2 cacheId, int2 cacheBlockOffset)
{
    // 2x2テクセル範囲を線形サンプル
    float2 fetchUV = (float2(cacheId*2 + cacheBlockOffset) + 1.0f) / TEXTURE_SIZE;
    fetchUV = clamp(fetchUV, 0.0f, 1.0f);
    s_aDepthCache[cacheId.x + cacheId.y*BLOCK_SIZE*2] = SampleLod0(g_DepthTexture, g_LinearSampler, fetchUV).r;
}

/// @brief  共有メモリからサンプルUV座標の位置に対応するキャッシュを取得
float GetDepthCache(float2 uv, int2 cacheBlockOffset)
{
    // 2x2テクセルを1つの値としてキャッシュしていることを考慮してキャッシュIDを計算
    uint2 cacheId = uint2(floor(uv * TEXTURE_SIZE - cacheBlockOffset));
    cacheId = clamp(cacheId/2, 0, BLOCK_SIZE * 2 - 1);
    return s_aDepthCache[cacheId.x + cacheId.y*BLOCK_SIZE * 2];
}

#endif // #ifdef USE_COMPUTE_SHADER

/// @brief AO計算時のView空間Position取得
float3 GetViewPosition(float2 uv, int2 cacheBlockOffset)
{
#ifdef USE_COMPUTE_SHADER
    const float depth = GetDepthCache(uv, cacheBlockOffset);
#endif // #ifdef USE_COMPUTE_SHADER
#ifdef USE_PIXEL_SHADER
    const float depth = SampleLod0(g_DepthTexture, g_LinearSampler, uv).r;
#endif // #ifdef USE_PIXEL_SHADER
    return DepthToViewPosition(depth, TexCoordToScreenSpace(uv), DEPTH_FACTOR_1, DEPTH_FACTOR_2, DEPTH_FACTOR_3, DEPTH_FACTOR_4);
}

/// @brief ワールド空間半球内部のランダムな点をサンプルし,XY空間に射影
float2 GetRamdomHemisphereXY(float2 uv, int seed)
{
    const float r1 = rand(uv, seed);
    const float r2 = frac(r1 * 2558.54f);
    const float r3 = frac(r2 * 1924.19f);
    const float r = pow(r1, 1.0f / 3.0f) * sqrt(1.0f - r3*r3);
    const float theta = 2 * PI * r2;
    float s, c;
    sincos(theta, s, c);
    return float2(c,s) * r;
}

/// @brief AOの計算
float CalculateAlchemyAO(float2 uv, int2 cacheBlockOffset = int2(0, 0))
{
    const float3 normal = DecodeNormal(SampleLod0(g_NormalTexture, g_LinearSampler, uv).rgb, NORMFMT_UNORM); // FIXME
    const float3 position = GetViewPosition(uv, cacheBlockOffset);

    /// ワールド空間でのサンプリング半径から,スクリーンスペース半径への変換係数
    // DEPTH_FACTOR_4 = 1/tan(fov/2)
    const float ry = DEPTH_FACTOR_4 * SAMPLE_RADIUS / position.z; // 画面に垂直な円のY軸半径
    const float2 radiusScale = ry * 0.5f * float2(ASPECT_RATIO, 1.0f); // 楕円スケール,0.5倍はUV空間へのマップ

    // AO値の積分
    float sum = 0.0f;
    for (int i = 0; i < SAMPLE_NUM; ++i)
    {
        const float2 sampleUV = uv + GetRamdomHemisphereXY(uv, i) * radiusScale; // サンプルするUV座標値を計算
        const float3 samplePosition = GetViewPosition(sampleUV, cacheBlockOffset); // view空間でのサンプリング地点
        const float3 v = samplePosition - position;
        sum += max(0.0f, dot(v, normal) + position.z*BIAS) / (dot(v, v) + EPSILON);
    }
    float AO = 2.0f * INTENSITY * sum / SAMPLE_NUM;

    // 遠くのAOは薄くするための,視錐台内部のz正規化値,これがないと無限遠方がおかしくなる
    const float depthAttenuate = saturate((FAR_Z - position.z) / (FAR_Z - NEAR_Z));
    AO *= depthAttenuate;

    return 1.0f - AO;
}

#ifdef USE_PIXEL_SHADER
float4 PS_Create(PPFX_OUT In) : SV_TARGET0
{
    return CalculateAlchemyAO(In.UV.xy);
}
#endif // #ifdef USE_PIXEL_SHADER

#ifdef USE_COMPUTE_SHADER
/// @brief エントリ
[numthreads(BLOCK_SIZE, BLOCK_SIZE, 1)]
void CS_Create(uint3 dispachThreadID : SV_DispatchThreadID, uint3 groupThreadID : SV_GroupThreadID, uint3 groupID : SV_GroupID )
{
    // キャッシュブロックのオフセット
    // 2x2テクセルごとに線形サンプルして,半径64pixel分のサンプルまで稼ぐ
    const int2 cacheBlockOffset = int2(groupID.xy) * BLOCK_SIZE - BLOCK_SIZE*3/2;

    // ViewPositionのキャッシュ読み込み
    LoadDepthCache(groupThreadID.xy, cacheBlockOffset);
    LoadDepthCache(groupThreadID.xy + int2(BLOCK_SIZE, 0), cacheBlockOffset);
    LoadDepthCache(groupThreadID.xy + int2(0, BLOCK_SIZE), cacheBlockOffset);
    LoadDepthCache(groupThreadID.xy + int2(BLOCK_SIZE, BLOCK_SIZE), cacheBlockOffset);
    GroupMemoryBarrierWithGroupSync();
    if ((dispachThreadID.x >= TEXTURE_SIZE_X) || (dispachThreadID.y >= TEXTURE_SIZE_Y)) return;

    // AO値の取得
    const float2 uv = (float2(dispachThreadID.xy)+ 0.5f) / TEXTURE_SIZE;
    g_DstTexture[dispachThreadID.xy].x = CalculateAlchemyAO(uv, cacheBlockOffset);
}
#endif // #ifdef USE_COMPUTE_SHADER

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

はじめに

前回,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

FFTの実装

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

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

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

クリックで展開

// #define INVERSE      :  逆フーリエ変換
// #define PASS_1ST        :  1pass目, X方向変換
// #define PASS_2ND        :  2pass目, Y方向変換

#if defined(PASS_1ST) && !defined(INVERSE)
#define INPUT_REALPART     // 入力は実数空間のみ
#endif // #if defined(PASS_1ST) && !defined(INVERSE)
#if defined(PASS_2ND) && defined(INVERSE)
#define OUTPUT_REALPART    // 出力は実数空間のみ
#endif // #if defined(PASS_1ST) && !defined(INVERSE)

// 共有メモリが許すならば,ダブルバッファを行う
#if LENGTH<=512
#define DOUBLE_BUFFER
#endif // #ifLENGTH<=512

Texture2D<float3> g_SrcTextureR   : register(t0);
#if !defined(INPUT_REALPART)
Texture2D<float3> g_SrcTextureI   : register(t1);
#endif // #if !defined(INPUT_REALPART)
RWTexture2D<float3> g_DstTextureR : register(u0);
#if !defined(OUTPUT_REALPART)
RWTexture2D<float3> g_DstTextureI : register(u1);
#endif // #if !defined(OUTPUT_REALPART)

static const uint REAL_PART = 0;
static const uint IMAGINARY_PART = 1;

///< @brief バタフライ演算用の計算バッファ
#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

/// @brief バタフライ演算に必要な重みと加算インデックスを取得
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);
    }
}

/// @brief バタフライ演算
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 // #ifndef DOUBLE_BUFFER

    // 加算
#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
/// @brief 逆フーリエ変換出力用のバタフライ演算パス
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 // OUTPUT_REALPART


/// @brief CSエントリポイント
[numthreads(LENGTH, 1, 1)]
void CS_Transform(uint3 position : SV_DispatchThreadID)
{
    const uint bufferId = position.x;
    // フェッチするテクスチャ座標の取得
#ifdef PASS_1ST
    // X方向
    uint2 texturePos = position.xy;
#endif // #ifdef PASS_1ST
#ifdef PASS_2ND
    // Y方向
    uint2 texturePos = position.yx;
#endif // #ifdef PASS_2ND
    texturePos = (texturePos + LENGTH / 2) % LENGTH; // 中心を周波数空間原点とする処理

    // 入力テクスチャを格納
    TempBuffer(0, bufferId, REAL_PART) = g_SrcTextureR[texturePos];
#if defined(INPUT_REALPART)
    TempBuffer(0, bufferId, IMAGINARY_PART) = 0.0f; // FFT入力は実部のみ
#else 
    TempBuffer(0, bufferId, IMAGINARY_PART) = g_SrcTextureI[texturePos];
#endif // #if defined(INPUT_REALPART)

    // バタフライ演算
    // 最後の1回のパスの計算は行われない,(直接RWTextureに書き出すため)
    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));
    }

    // バタフライ演算最後の1回
    // 出力へ書き出し
    {
        GroupMemoryBarrierWithGroupSync();
        const uint srcTempId = (BUTTERFLY_COUNT - 1)% 2;
#ifdef OUTPUT_REALPART
        InverseButterflyPassRealPart(BUTTERFLY_COUNT - 1, bufferId, srcTempId, g_DstTextureR[texturePos]); // IFFT出力は実部のみ
#else
        ButterflyPass(BUTTERFLY_COUNT - 1, bufferId, srcTempId, g_DstTextureR[texturePos], g_DstTextureI[texturePos]);
#endif // #ifdef OUTPUT_REALPART
    }
}

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

FFTの適用結果

512x512サイズの画像に対して,FFTを適用しました。 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ができそうです。

Compute ShaderでSummed Area Table

Summed Area Table

コンピュートシェーダーに興味が出てきたので,練習がてらSummed Area Table作成してブラーをかけてみます。

https://docs.nvidia.com/gameworks/content/gameworkslibrary/graphicssamples/d3d_samples/d3dcomputefiltersample.htm

https://www.slideshare.net/egorodet/cpu-is-in-focus-again-implementing-dof-on-cpu

Summed Area Tableとは,テクスチャのテクセルと表したときに,以下のようにX方向,Y方向に順番に和を取っていったテーブルのことです

Summed Area Tableを作成しておくと,テーブルから4点フェッチして重みを工夫してサンプルするだけで, 元のテクスチャの任意の矩形範囲内部のピクセルの積分が計算できます。

テーブルさえ作ってしまえば,どんな積分範囲であっても高速に計算ができるので, 被写界深度処理のような ピクセルによってブラー範囲(=積分範囲)を変えるというケースで特に有効かもしれません。

ただし,Summed Area Table自体の作成するときに, 全ピクセルに対して依存するピクセルの総和と計算しておく必要があるので,この計算コストが問題となります。

Summed Area Tableの作成

Summed Area Tableの作成には,Compute Shaderを利用しました。
実際のシェーダーコードとCPU側のソースは次のような感じです。

//// HLSL側のCompute Shader
Texture2D<float4> g_SrcTexture : register(t0);
SamplerState g_LinearSampler : register(s0);
RWTexture2D<float4> g_DstTexture : register(u0);
#define    PIXEL_SIZE g_aTemp[0].xy
groupshared float3 s_aTempBuffer[LENGTH];

[numthreads(LENGTH, 1, 1)]
void CS_CreateTable(uint3 position : SV_DispatchThreadID)
{
    const uint bufferId = position.x;
#ifdef PASS_1ST
    // 水平パス
    const uint2 texturePos = position.xy;
#endif
#ifdef PASS_2ND
    // 垂直パス
    const uint2 texturePos = position.yx;
#endif
        // テクスチャの読み込み
    s_aTempBuffer[bufferId] = SampleLod0(g_SrcTexture, g_LinearSampler, float2(texturePos)*PIXEL_SIZE).rgb;

    // イテレーション
    float3 addBuf;
    for (uint i = 0; i < ITERATION; ++i)
    {
        GroupMemoryBarrierWithGroupSync();
        const uint delta = (0x1 << i);
        const float3 buffer = (bufferId >= delta) ? (s_aTempBuffer[bufferId - delta]) : (0.0f);
        GroupMemoryBarrierWithGroupSync();
        s_aTempBuffer[bufferId] += buffer;

    }
    float4 ret = float4(0.0f, 0.0f, 0.0f, 1.0f);
    ret.rgb = s_aTempBuffer[bufferId];
    g_DstTexture[texturePos] = ret;
}

//// CPP側の擬似コード
// 水平方向加算
auto pPass = rPfxEffect.GetPass("SummedAreaTableCreate", "L512_1ST").lock();
auto pTempBuffer = rPfxEffect.GetTempTexturePool().Get(dstSize, dstSize, tableFormat, bindFlags).GetTextureView();
pPass->CSBindConstantBuffer(0, tempCB);
pPass->CSBindTexture(0, pSrcTextureView);
pPass->CSBindSampler(0, rPfxEffect.GetLinearSampler());
pPass->CSBindRWTexture(0, pTempBuffer);
result = rPfxEffect.Dispatch(pDeviceContext, pPass, 1, dstSize, 1);
// 垂直方向加算
pPass = rPfxEffect.GetPass("SummedAreaTableCreate", "L512_2ND").lock();
auto pResultTextureView = rPfxEffect.GetTempTexturePool().Get(dstSize, dstSize, tableFormat, bindFlags).GetTextureView();
pPass->CSBindConstantBuffer(0, tempCB);
pPass->CSBindTexture(0, pTempBuffer);
pPass->CSBindSampler(0, rPfxEffect.GetLinearSampler());
pPass->CSBindRWTexture(0, pResultTextureView);
result = rPfxEffect.Dispatch(pDeviceContext, pPass, 1, dstSize, 1);

このコードの要点をまとめると

• X方向の加算とY方向の加算を別パスに分離
• 再帰的に加算をして計算量をテクスチャ幅のlogに落とす
• シェーダーのグループ共有メモリを利用して,テクスチャフェッチ回数を最小限に減らす

という感じです。

何も考えずに加算をすると,テクスチャ幅]として,一つのピクセルあたり回のループを回すことになります。
これを,以下の図のように2べき飛ばしで再帰的に加算をしていくことで, 回まで減らすことが可能です。

この再帰計算を行うときに便利なのがgroupsharedで宣言されたグループ共有メモリです。 個のスレッドを走らせるグループを個作成し,それぞれのグループで個のfloat3を計算バッファとして確保しておきます。
最初の一回だけこのバッファに対してテクスチャ読み出しをおこなって,あとの再帰的な加算を行うときはバッファを利用することで,テクスチャ読み出しを必要最低限に抑えます。

Summed Area Tableの作成自体はVertex ShaderとPixel Shaderでも工夫をすれば可能ですが,再帰的に加算をする度にテクスチャ読み出しをするので, Compute Shaderほど高速には行なえません。 DirectX11ではグループ共有メモリは32KBまでしか利用できないので,float4で考えるとテクスチャサイズは2048が限界ということに注意が必要です。 また,加算時の誤差が激しいので,32bit浮動小数テクスチャを利用しています。

ブラー結果

512x512ピクセルのSummed Area Tableを作成すると,次のような感じになりました。

入力テクスチャ	Summed Area Table

Summed Area Tableの右端の値は9000ぐらいで,普通に表示すると飽和して見えないので,適当にスケールしてあります。
このテーブルを利用して,ただ矩形状に平均値を計算するボックスフィルタを適用してみます。

だいぶカクカクとした結果になりました。 今度は矩形を3つ重ねて少し滑らかになるようなフィルタを適用してみます。

まずまずではないでしょうか。
このパス全体のGPU実行時間は,GTX970の環境で

• 512x512テクスチャのTable作成 : 0.14 ms
• 1280x720のテクスチャへのブラー書き出し : 0.34 ms

程度のGPU時間でした。なかなか高速だと思います。

Log Summed Area Table

おまけです。
Summed Area Tableを作成するときに,対数平均を取るようにして,ブラー時にexpで元に戻すようにしたら,誤差が減るのではないかと思ってやってみました。

Summed Area Table	Log Summed Area Table

Summed Area Table	Log Summed Area Table

少しブラーの角が取れたような気もします。