アナリティック・フォグレンダリング：ボリュメトリックプリミティブを用いた手法（2025）

この投稿では、密度が変化する霧の領域を描画するためのシンプルなテクニックについて解説します。
まずは霧レンダリングの基本原理とよく使われる手法をざっくり紹介し、その後に具体的なテクニックを説明します。数式が登場しますので、予めご注意ください！

---

1. 霧レンダリングの簡単レビュー

霧の描画は「光が媒質を通過する際にどのように振る舞うか」をシミュレートすることに帰着します。最も基本的で、今回焦点を当てるのは吸収（光エネルギーの損失）です。言い換えれば、光が通過する「もの」が多ければ多いほど、エネルギーは減衰します。

光線が媒質を通過するときに失うエネルギー

光線が媒質（影付け領域）を通過した後の出射光と入射光の比率を 透過率 (transmittance) と呼び、次式で定義されます。

[ T = \frac{L_{\text{out}}}{L_{\text{in}}} ]

ここで (L_{\text{in}}) は媒質に入る光量、(L_{\text{out}}) は出る光量です。
媒質の密度と透過率との関係はビール・ラマート定理（Beer–Lambert Law）で表されます。

[ T = e^{-\int_{0}^{s} \rho(\mathbf{x}(t)),dt} ]

(\rho) は光線上を通る総密度です。(T) がわかれば、単に (L_{\text{in}}) に掛けて出射光量が決まります。

媒質の密度が一定 (\rho_0) であれば、

[ T = e^{-\rho_0,s} ]

ここで (s) は光線内にある長さです。
この式はゲームなどでよく使われる「距離霧（distance fog）」の基礎です。時にはビール・ラマート定理を完全に省略し、単に総密度だけで遠距離のオブジェクトをフェードアウトさせます。いずれにしても、常に一定密度の媒質が全視域を覆っていると仮定しています。

簡易的な指数型距離霧の例

---

1.1 限界付き領域での柔軟性

媒質を境界体（plane, sphere, box 等）内に制限すれば、霧の見え方をより細かく調整できます。視線がこれらのプリミティブと交差する点を求め、その内部長さ (s) を計算し、同じ式で総密度を求めます（ここでも一定密度を仮定）。

球・箱・平面を使って霧領域を囲む方法

（光線―プリミティブの交差判定は本投稿では省略。詳細は Inigo Quilez のウェブサイト参照）

---

2. 異種体（heterogeneous volumes）

実際、雲や煙などの媒質は全域で一定密度になることは稀です。密度が変化している場合を 異種体 と呼びます。このような媒質では、より多様なエフェクトが可能ですが、式も複雑になります。

透過率は単なる積ではなく、光線上での積分となります。

[ T = e^{-\int_{s_0}^{s_1} \rho(\mathbf{x}(t)),dt} ]

ここで (s_0, s_1) は境界体内にある光線の始点・終点です。関数 (\rho(\cdot)) は空間上の任意点での密度を表します。

多くの場合、閉形式積分が得られないため、数値的手法（小ステップで密度をサンプリングし、矩形近似で合計）に頼ります。これが実際に行っている レイマーチング です。非常に柔軟で、任意の位置から (\rho) をサンプリングできれば（例：ボリュームテクスチャ）、積分を数値的に評価できます。

ただし、サンプル数が少ないとエイリアシングが発生します。解析解は柔軟性が低いものの、エイリアシングがなく高速な場合があります。そのため、異なる媒質ごとの解析式を求めることに意義があります。

---

3. 放射関数（radial functions）

ある点 (\mathbf{p}) が原点からどれだけ離れているかにのみ依存する関数は 放射関数 と呼ばれます。これらは扱いやすく、密度が領域内で変化するボリュームを定義できます。複数の放射密度関数を組み合わせることで、より複雑な霧や煙エフェクトを作り出せます。

単純な放射関数の場合、線分上での積分はそれほど難しくありません。放射関数は距離 (r=|\mathbf{p}|) に依存します。光線が原点 (\mathbf{o})、方向ベクトル (\mathbf{d}) を持ち、距離 (t) である点は

[ \mathbf{p}(t)=\mathbf{o}+t,\mathbf{d} ]

よって光線上の任意点への距離は

[ r(t)=|\mathbf{o}+t,\mathbf{d}| ]

与えられた放射関数を (t) の関数に書き換えると、新しい関数が得られ、これを積分できれば光線上の密度変化が分かります。これは元のボリュームの「横断面」を表すものです。

光線方程式に線形変換を掛けることで、任意にスケール・回転したボリュームを配置できます（変換行列 (M) を用いて）

[ \mathbf{p}(t)=M,(\mathbf{o}+t,\mathbf{d})=\mathbf{o}'+t,\mathbf{d}' ]

残りは新しい関数の不定積分を求め、光線始点・終点に評価するだけです。
もし (\mathbf{o}) がプリミティブから遠い場合は、数値精度問題が生じることがあります。その際は (t=0) を区間の開始点とし、新しい限界を 0 と (\Delta t = t_1-t_0) に置き換えることで安定化します。

---

4. 放射関数例と不定積分

以下に、いくつかの放射密度関数、その光線上での横断面、および解析的不定積分を示します（GLSL/GLSL‑style のコードスニペットを想定）。

4.1 線形密度 (Triangular in 1D / cone in 2D)

[ \rho(r)=\max(0,;1-a r^2 + b r + c) ]

横断面:

[ f(t)=\max(0,;1-a t^2 + b t + c) ]

不定積分（平方完成後）：

float antiderivativeLinear(float a, float b, float c, float t)
{
    float u = t + b / a;
    float v = (b * b - a * c) / (a * a);
    float radical = sqrt(u * u + v);

    return t - sqrt(a) * (u * radical - v * log(0.001f + u + radical)) / 2.0f;
}

4.2 二次密度 (Parabolic shape)

[ \rho(r)=\max(0,;1-a r^4 + b r^3 + c r^2) ]

横断面:

[ f(t)=\max(0,;1-a t^4 + b t^3 + c t^2) ]

不定積分（多項式積分）：

float antiderivativeQuadratic(float a, float b, float c, float t)
{
    return t * (t * (t * a / -3.0f - b) - c + 1.0f);
}

4.3 四次密度 (Smoothstep‑like shape)

[ \rho(r)=\max(0,;1-a r^6 + b r^5 + c r^4) ]

横断面:

[ f(t)=\max(0,;1-a t^6 + b t^5 + c t^4) ]

不定積分：

float antiderivativeQuartic(float a, float b, float c, float t)
{
    return 
        t * (t * (t * (t * (a * a * t / 5.0f + a * b) + 
        (2.0f * b * b + a * c - a) * 2.0f / 3.0f) +
        (b * c - b) * 2.0f) + (c * c - 2.0f * c + 1.0f));
}

4.4 スパイキー密度 (Peaks sharply in the centre)

[ \rho(r)=\max(0,;1-a r^8 + b r^7 + c r^6) ]

横断面:

[ f(t)=\max(0,;1-a t^8 + b t^7 + c t^6) ]

不定積分（展開・簡略化後）：

float antiderivativeSpiky(float a, float b, float c, float t)
{
    float u = t + b / a;
    float v = (b * b - a * c) / (a * a);
    float radical = sqrt(u * u - v);

    return 
        t * (t * (t * a / 3.0f + b) + c + 1.0f) - 
        sqrt(a) * (u * radical - v * log(0.001f + u + radical));
}

4.5 ガウス密度 (Classic bell‑curve)

[ \rho(r)=\exp(-a r^2 + b r + c) ]

横断面:

[ f(t)=\exp(-a t^2 + b t + c) ]

不定積分は誤差関数 (\operatorname{erf}) を用います：

float erf(float z)
{
    float az = abs(z);
    return sign(z) * (1.0f - 1.0f / pow(1.0f + az * (0.278393f + az *
        (0.230389f + az * (0.000972f + az * 0.078108f))), 4.0f));
}

float antiderivativeGaussian(float a, float b, float c, float t)
{
    float sqrtA = sqrt(a);
    return sqrt(PI) * exp(b * b / a - c) *
           erf((a * t + b) / sqrtA) / (2.0f * sqrtA);
}

---

5. 補足事項

5.1 ライティングとシャドウ

これらのボリュームが光や影にリアルに反応するようにするには、完全な物理方程式は解析解で求めるのが難しいため、各プリミティブに簡易シェーディングモデルを適用するというハックがあります。全体として媒質が期待通りに振舞うように見えます。

5.2 動的霧とブーリアン演算

ボリュームをプリミティブで表現すると、さまざまな応用が可能です。たとえば、パーティクルシステムと連携して動的・アニメーション化された霧を作ることもできますし、布尔演算でより複雑な境界体（例：Fog of War）を構築することもできます。

---

このテクニックが興味深いと思われたら、ぜひプロジェクトに取り入れてみてください。全関数の可視化とソースコードは私の Shadertoy 実装で確認できます。ご覧いただきありがとうございました！