流体シミュレーションの最適化と描画技術

はじめに

こんにちは。株式会社スパーククリエイティブCTOの広本です。
普段は SPARK GEARの開発をしたり、UEでエンジンの改造をしたり、Unityでレンダラーや揺れもの物理を作ったりしています。

先日「0から始める流体シミュレーション」という内容で社内勉強会に参加しました。
80ページ以上の資料を作ったのですがせっかくなのでブログにも残しておこうと思います。

今回作成した流体シミュレーションです。

流体シミュレーションにおいて、計算そのものの精度（ナビエ・ストークス方程式の解法など）はもちろん重要ですが、最終的なアウトプットの質を左右するのはやはり描画技術です。

本記事では、ボリューメトリックな流体を描画する際に直面した課題と、それをどのように解決していったのかを雑に解説していきます。

はじめに
1. 「レイマーチング」の実装
2. 「ジッタリング」の実装
3. 「TAA」の実装
4. 「ゴースティング」の問題
5. 「遮蔽ボリューム」によるセルフシャドウ
- 遮蔽ボリュームの事前計算
- ジッタリングとブラーによる「ソフトシャドウ」化
6. 「散乱」の表現
7. Cubicフィルタリング
8. まとめ
おまけ

1. 「レイマーチング」の実装

まずは、流体ボリュームを描画するための基礎技術としてレイマーチングを実装しました。

レイマーチングは、カメラから視線方向（レイ）を飛ばし、一定の間隔でボリュームデータをサンプリングして密度を積算していく手法です。
一般的なポリゴン描画とは異なり、煙や炎といった「形が不定で内部に密度を持つ物体」を描画するのに適しています。

というわけで、シンプルなレイマーチングの実装の結果がこちらです。

見ての通り激しい「モアレ（縞模様）」が発生してしまいました。
これは、レイを進めるステップ間隔が固定であるために発生するサンプリングアーティファクトです。

ステップ数を極端に増やせばある程度は解消できますが、パフォーマンスが死ぬので出来ればステップ数を増やさずに改善したいところです。

2. 「ジッタリング」の実装

モアレが発生する根本的な原因は、サンプリング位置が規則正しすぎることにあります。

そこで導入したのがジッタリングです。

これは、レイの開始位置をピクセルごとにランダムにずらすという手法です。

これにより、規則的な縞模様（エイリアシング）を、不規則な「ノイズ」へと変換することができます。

コード抜粋

float Jitter = hash(ScreenUV + cbJitterOffset);
RayStart += Direction * Jitter * cbJitteringStrength;

というわけで、ジッタリングの実装するとこんな感じになります。
ここでは適当なノイズを使っていますが、実際の実装ではブルーノイズを使ってジッタリングを行っています。

一応不快なモアレは綺麗に消えました。

しかし、その代償として画面全体が非常にノイジーになってしまいました。

モアレに比べればノイズの方がマシだと思いますがこれだけでは問題の種類が変わっただけに過ぎません。

3. 「TAA」の実装

ジッタリングによって生じた高周波ノイズを取り除くために採用したのが、TAA（Temporal Anti-Aliasing）です。

TAAは、現在のフレームの情報だけでなく、過去のフレームの結果をうまくブレンドすることで、時間軸方向でサンプリング数を稼ぐ手法です。

ジッタリングで毎フレーム異なる位置をサンプリングしているため、時間をかけてそれらを平均化すれば、あたかも「超高解像度でサンプリングした」かのような滑らかな結果が得られるはずです。

コード抜粋

float4 main(VS_PS_POSTFX input) : SV_Target
{
    float2 uv = input.Texcoord0;

    // 現在のフレームの色を取得
    float4 CurrentColor = g_Texture0.SampleLevel(g_Sampler0, uv, 0);

    // カメラ行列を使って前のフレームのUV座標を計算
    // a. 簡易版: 現在のクリップ座標を使う
    float4 CurrentClipPos = float4(uv * float2(+2.0, -2.0) + float2(-1.0, +1.0), 0.0, 1.0);

    // b. 前のフレームでのクリップ座標を計算
    float4 PrevClipPos = mul(CurrentToPrevViewProj, CurrentClipPos);

    // c. クリップ座標を正規化デバイス座標に変換し、UV座標に変換
    PrevClipPos.xyz /= PrevClipPos.w;
    float2 prev_uv = PrevClipPos.xy * float2(+0.5, -0.5) + 0.5.xx;

    // ヒストリをサンプリング
    float4 HistoryColor;
    if ((prev_uv.x >= 0.0) && (prev_uv.x <= 1.0) && (prev_uv.y >= 0.0) && (prev_uv.y <= 1.0))
    {
        HistoryColor = g_Texture1.SampleLevel(g_Sampler1, prev_uv, 0);
    }
    else
    {
        HistoryColor = CurrentColor; // 画面外なら現在の色を使う (ゴースト対策)
    }

    // 現在のピクセルの3x3近傍から色のAABB(最小値/最大値)を計算
    float4 MinColor = CurrentColor;
    float4 MaxColor = CurrentColor;
    float2 TexelSize = 1.0.xx / ViewportSize.xy;
    [unroll]
    for (int y = -1; y <= +1; ++y)
    {
        [unroll]
        for (int x = -1; x <= +1; ++x)
        {
            float4 NeighborColor = g_Texture0.SampleLevel(g_Sampler0, uv + float2(x, y) * TexelSize, 0);
            MinColor = min(MinColor, NeighborColor);
            MaxColor = max(MaxColor, NeighborColor);
        }
    }

    // ヒストリの色をクランプ
    HistoryColor = clamp(HistoryColor, MinColor, MaxColor);

    // ブレンド
    float4 FinalColor = lerp(HistoryColor, CurrentColor, 0.1);

    return FinalColor;
}

UE界隈の方々にはおなじみですが実装結果はこんな感じです。

ノイズが消え、非常に滑らかで高品質な表現が可能になりました。

レイマーチングの重厚感を保ちつつ、リアルタイムでここまでの品質が出せればひとまずOKでしょう。

今回はメッシュの描画にもAAがかかるようにハルトンシーケンスを利用したカメラ位置のジッタリングも行っています。

コード抜粋

// ハルトンシーケンスの1次元の値を計算する関数
float Halton(int index, int base)
{
    float result = 0.0f;
    float f = 1.0f / (float)base; // 基数の逆数 (1/b)

    while (index > 0)
    {
        result += (float)(index % base) * f; // 余り * (1 / b^k) を加算
        index /= base;                       // indexを基数で割る
        f /= (float)base;                    // 次の桁の重み (1 / b^(k+1))
    }

    return result; // 0.0 から 1.0 の間の値
}

しかし、TAAの導入により、新たな、そして根深い課題が浮上することになります。

4. 「ゴースティング」の問題

TAAは「過去のフレーム」を利用する技術です。

そのため、画面内で物体が激しく動いた場合、現在の位置と過去の位置にズレが生じます。

通常はモーションベクトルを使って補正を行いますが、流体シミュレーション、特に半透明な煙や炎の場合、これが非常に厄介です。

こちらもUE界隈の方々にはおなじみのゴースティング（残像現象）です。

特に動きの速い部分において、過去のフレームの絵が亡霊のように残ってしまっています。

半透明描画は深度（Depth）があいまいであったり、流体のように形状そのものが毎フレーム変形するオブジェクトの場合、正確な履歴の追跡が困難であるため、TAAとの相性は本質的に悪いです。

ここに関しては効果的な解決方法も無いのとカメラを高速で動かすケースはあまり想定していないので諦めました。

5. 「遮蔽ボリューム」によるセルフシャドウ

TAAによって滑らかな描画は手に入りましたが、この段階ではまだ煙が「のっぺり」としており、立体感が不足しています。

これを解決するのが、光の濃淡計算である「セルフシャドウ」の実装です。

遮蔽ボリュームの事前計算

セルフシャドウの実装として、レイマーチングの各ステップからさらに光源方向へレイを飛ばす手法もありますが、これは計算負荷が爆発的に増大してしまいます。

そこで今回は、事前に「遮蔽マップ」をボリュームデータとして計算しておく手法を採用しました。

遮蔽情報の生成: メインの描画パスの前に、光源方向からの光の透過具合を計算し、別の3Dテクスチャ（ボリューム）に焼き込みます。
描画時の参照: メインのレイマーチング時は、その3Dテクスチャをサンプリングするだけで、「その地点に光がどれだけ届いているか」を高速に取得できます。

というわけで実装してみました。


▲ セルフシャドウなし	▲ セルフシャドウあり

ジッタリングとブラーによる「ソフトシャドウ」化

しかし、単純に遮蔽ボリュームを計算しただけでは、影のエッジがクッキリしすぎてしまい、煙らしい柔らかさが損なわれることがあります。

そこで、より自然な「ソフトシャドウ」を実現するために、ここでもジッタリングとブラーのテクニックを導入しました。

光源方向のジッタリング
遮蔽ボリュームを生成する際、光源に向かうレイの方向をピクセルごとに微小にランダムにずらします。

これにより、影のエッジが適度に散乱し、擬似的なソフトシャドウ効果が得られます。

ただし、このままでは影がノイズだらけになってしまいます。

遮蔽ボリュームへのガウスブラー
生成されたノイズ混じりの遮蔽ボリュームに対して、ガウスブラーフィルタを適用します。

コード抜粋

[numthreads(8, 8, 8)]
void cs_main(uint3 thread_id : SV_DispatchThreadID)
{
    int3 index = (int3)thread_id;
    float3 uvw = mad(index, cbRecipShadowDimensionSize.xyz, cbHalfShadowVoxelSize.xyz);

    const int BlurRange = 1;
    const float uvwRange = cbShadowBlurRange;
    float totalWeight = 0;
    float LightPower = 0;
    for (int x = -BlurRange; x <= BlurRange; x++) {
        for (int y = -BlurRange; y <= BlurRange; y++) {
            for (int z = -BlurRange; z <= BlurRange; z++) {
                float3 offset = float3(x, y, z) * uvwRange.xxx;
                float3 samplePos = uvw + offset;
                samplePos = saturate(samplePos);
                float dist = length(offset);
                float w = gaussian(dist);
                totalWeight += w;
                LightPower = mad(g_TextureShadow.SampleLevel(g_SamplerShadow, samplePos, 0.0).r, w, LightPower);
            }
        }
    }
    LightPower /= totalWeight;
    g_TextureShadowRW[index] = LightPower;
}

これにより、ノイズが滑らかに平均化され、ソフトシャドウが生成されます。


▲ ソフトシャドウなし	▲ ソフトシャドウあり

この「遮蔽ボリューム」と「ソフトシャドウ化」のプロセスを経ることで、煙の奥まった部分に落ちる影が柔らかく広がり、リアルタイム性を損なうことなく、レンダリング時にソフトシャドウを実現できました。

6. 「散乱」の表現

影を落とすだけでは、単に黒い塊になってしまいます。

現実の煙や雲は、光を受けると内部で光が拡散し、柔らかく輝いて見えます。

これを再現するのが散乱（Scattering）の計算です。

物理的に厳密な計算（多重散乱など）を行うと重すぎるため、簡易的なライティングモデルを採用しています。

Beer-Lambertの法則: 密度に応じた光の減衰
Henyey-Greenstein位相関数: 光が入射角に対してどの方向に散乱するか（逆光で縁が光る表現など）

コード抜粋

float HenyeyGreensteinPhaseFunction(float cosTheta, float g)
{
    static const float INV_FOUR_PI = 1.0 / (4.0 * 3.14159265);
    float g2 = g * g;
    float numerator = 1.0 - g2;
    float denominator = pow(max(0.0, 1.0 + g2 - 2.0 * g * cosTheta), 1.5);
    // ゼロ除算を避ける
    denominator = max(denominator, 1e-6);
    return INV_FOUR_PI * numerator / denominator;
}

というわけで実装前後の比較画像です。


▲ 散乱なし	▲ 散乱あり

これらを組み合わせることで、ただの煙の塊が「光を透過・拡散する物質」として認識されるようになります。

7. Cubicフィルタリング

最後に、地味ながらも映像品質に直結する「フィルタリング」の問題について触れておきます。

流体シミュレーションは計算コストが高いため、描画解像度（画面サイズ）よりも低い解像度のグリッド（例：128x128x128など）で計算を行うのが一般的です。

しかし、これを単純に拡大して描画（Linear補間）すると、グリッドのマス目が目立ってしまい、カクカクした見た目になってしまいます。

そこで、テクスチャサンプリング時にCubicフィルタリングを導入しました。

コード抜粋

float2 GetReactionYBSlpline(float3 uvw, float OffsetCoef, float3 HG_z) {
    float3 uvw_y = uvw + float3(0, OffsetCoef, 0) * cbRecipDimensionSize.y;

    float3 uvw_z0 = uvw_y + float3(0, 0, HG_z.x) * cbRecipDimensionSize.z;
    float3 uvw_z1 = uvw_y - float3(0, 0, HG_z.y) * cbRecipDimensionSize.z;

    float2 Reaction0 = g_TextureReaction.SampleLevel(g_SamplerReaction, uvw_z0, 0.0);
    float2 Reaction1 = g_TextureReaction.SampleLevel(g_SamplerReaction, uvw_z1, 0.0);
    return lerp(Reaction0, Reaction1, HG_z.z);
}

float2 GetReactionXBSlpline(float3 uvw, float OffsetCoef, float3 HG_y, float3 HG_z) {
    float3 uvw_x = uvw + float3(OffsetCoef, 0, 0) * cbRecipDimensionSize.x;
    float2 A = GetReactionYBSlpline(uvw_x, +HG_y.x, HG_z);
    float2 B = GetReactionYBSlpline(uvw_x, -HG_y.y, HG_z);
    return lerp(A, B, HG_y.z);
}