はじめに

この記事は球面調和関数とCG. Part1, ~Light Probeは何を計算しているのか？~の続きになります。前回解説した理論を参考にして、実際に skybox に対する light probe を実装してみようと思います。

Light probe

まずは、light probe はどんな値を持っているのかおさらいしていきます。

Light probe の目的は diffuse 反射の irradiance を記録することです。法線 $n$ に対する irradiance $E(n)$ は入射輝度 $L_i(\omega_i)$ に対して、

$$
E(n) = \int_{\Omega(n)} \frac{1}{\pi} (n \cdot \omega_i) L_i(\omega_i) d\omega_i
$$

と表されます（正確には irradiance を $1/\pi$ でスケールしたものです。ここではこれも irradiance と呼ぶことにします)。Light probe はこの球面上の関数 $E(n)$ を保持します。

$E(n)$ の持ち方は直接 texture として保持する irradiance environment map という方法もありますが、一般には SH による近似手法が使われています。

SH $Y_{\ell,m}$ は球面上の関数をフーリエ級数のように関数展開することができます。この性質を利用して、低次の SH によって irradiance $E(n)$ を近似します。

$$
E(n) \approx E_{0,0} Y_{0,0}(n) + E_{1,-1} Y_{1,-1}(n) + E_{1,0} Y_{1,0}(n) + …
$$

この時、$E(n)$ の計算に必要なパラメーターは各 SH $Y_{\ell,m}$ の係数 $E_{\ell,m}$ だけになります。したがって、irradiance environment map より少ないパラメーターで $E(n)$ を保持することが可能となります。これが light probe における SH 圧縮の仕組みです。

一般的には 2 次までの SH が採用されており(1 次までの場合も)、RGB それぞれで SH 係数を持つことから計 27 個のパラメーターが必要となります。これを事前計算します。

Irradiance における SH は畳み込みの性質から入射輝度にかかわる SH 係数 $L_{\ell,m}$ とそれ以外に関わる係数 $\hat{A}_{\ell}$ に分解することができます。

$$
E_{\ell,m} = \hat{A}_\ell L_{\ell,m}
$$

$$
L_{\ell,m} = \int_{S^2} L_i(\omega) Y_{\ell,m}(\omega) d\omega
$$

$\hat{A}_{\ell}$ は解析解があるため、私たちが事前計算すべきことは $L_{\ell,m}$ だけであり、周囲の輝度値さえわかれば数値計算することができます。

実装

今回は skybox に対する light probe を実装してみました。言語は C++ & Slang, APIは Vulkan を使用しています。また、今回の実装は理論との整合性のため、最適化はせずに分かりやすさを重視しています。

SH の数値積分は CPU 上で行いました。その疑似コードを以下に記載します。

std::array<float, 9> EvaluateSH9(const float3& d)
{
    const float x = d.x;
    const float y = d.y;
    const float z = d.z;

    constexpr float A0 = Pi<float>;
    constexpr float A1 = 2.0f * Pi<float> / 3.0f;
    constexpr float A2 = Pi<float> / 4.0f;

    const float C00 = A0 * 0.5f * sqrt(1.0f / Pi<float>);
    const float C1 = A1 * sqrt(3.0f / (4.0f * Pi<float>));
    const float C2 = A2 * 0.5f * sqrt(15.0f / Pi<float>);
    const float C20 = A2 * 0.25f * sqrt(5.0f / Pi<float>);
    const float C22 = A2 * 0.25f * sqrt(15.0f / Pi<float>);

    return {
        C00,
        C1 * y,
        C1 * z,
        C1 * x,
        C2 * x * y,
        C2 * y * z,
        C20 * (3.0f * z * z - 1.0f),
        C2 * x * z,
        C22 * (x * x - y * y),
    };
}
...

    float* pixel; // image data
    for (int i = 0; i < image.Width(); i++)
    {
        for (int j = 0; j < image.Height(); j++)
        {
            uint32_t index = image.Channel() * (i + image.Width() * j);
            float3 li = { pixel[index], pixel[index + 1], pixel[index + 2] };

            float u = (i + 0.5) / image.Width();
            float v = (j + 0.5) / image.Height();
            float theta = v * Pi<float>;
            float phi = u * Tau<float>;
            float3 d = { sin(theta) * cos(phi),
                         sin(theta) * sin(phi),
                         cos(theta) };

            const float dTheta = Pi<float> / image.Height();
            const float dPhi = Tau<float> / image.Width();
            const float solidAngle =
                std::sin(theta) * dTheta * dPhi;

            auto sh = EvaluateSH9(d);
            for (int i = 0; i < 9; i++)
            {
                shc[i] += li * sh[i] * solidAngle / Pi<float>;
            }
        }
    }

EvaluateSH9 は与えられた方向に対して、L0 から L2 までの SH を返す関数です。ここでは単純化するため、cosine 分布の係数である $A_\ell$ も含めて計算しています。

一方で、肝心の積分部分は下の方にある for 文の中で行われています。ここでは各 pixel に対して被積分関数と微小量の計算を行い、実際に数値積分を行っています。

まず $L_i$ ですが、これはシンプルに pixel の色です。コード上ではここに相当します

            float3 li = { pixel[index], pixel[index + 1], pixel[index + 2] };

次にやることは入射光である pixel の方向です。一般的な skybox に使われる HDRI は世界地図のように縦が緯度、横が経度となるように球面上にマッピングされています。なので、それに従うように緯度、経度を計算して入射方向の算出を行っています

            float u = (i + 0.5) / image.Width();
            float v = (j + 0.5) / image.Height();
            float theta = v * Pi<float>;
            float phi = u * Tau<float>;
            float3 d = { sin(theta) * cos(phi),
                         sin(theta) * sin(phi),
                         cos(theta) };

次は積分に使用する微小量の計算です。微小立体角 $d\omega$ は微小角 $d\theta, d\phi$ によって次のように計算されます。

$$
d\omega = \sin{\theta} d\theta d\phi
$$

この式に従い、微小立体角の計算をしているのがこの部分です。

            const float dTheta = Pi<float> / image.Height();
            const float dPhi = Tau<float> / image.Width();
            const float solidAngle =
                std::sin(theta) * dTheta * dPhi;

あとはこれらを合わせて積分値を計算すればOKというわけです。各 SH を取得し、求められた被積分関数と合わせて、最終的な結果を加算していきます。

            auto sh = EvaluateSH9(d);
            for (int i = 0; i < 9; i++)
            {
                shc[i] += li * sh[i] * solidAngle / Pi<float>;
            }

これでベイク側のやることは終わりました。次はこれを shader 側でデコードして irradiance を取得します。SH から irradiance を取得する最小限の shader の全文を以下に記載します。

[vk::binding(0)]
cbuffer Settings : register(b0)
{
    float3 shc00;
    float padding1;
    float3 shc1m1;
    float padding2;
    float3 shc10;
    float padding3;
    float3 shc1p1;
    float padding4;
    float3 shc2m2;
    float padding5;
    float3 shc2m1;
    float padding6;
    float3 shc20;
    float padding7;
    float3 shc2p1;
    float padding8;
    float3 shc2p2;
    float padding9;
};

struct PushConstant
{
    row_major float4x4 cameraMatrix;
};

[vk_push_constant]
PushConstant pushConstant;

struct VertexInput
{
    [[vk::location(0)]]
    float3 position : POSITION;
    [[vk::location(1)]]
    float3 normal : NORMAL;
};

struct VertexOutput
{
    float4 position : SV_Position;
    [[vk::location(0)]]
    float3 normal : TEXCOORD0;
};

static const float Pi = 3.14159265;
float3 EvaluateIrradianceSH(float3 d)
{
    d = normalize(d);

    const float x = d.x;
    const float y = d.y;
    const float z = d.z;

    const float C00 = 0.5f * sqrt(1.0f / Pi);
    const float C1 = sqrt(3.0f / (4.0f * Pi));
    const float C2 = 0.5f * sqrt(15.0f / Pi);
    const float C20 = 0.25f * sqrt(5.0f / Pi);
    const float C22 = 0.25f * sqrt(15.0f / Pi);

    float3 result = 0.0;
    result += shc00 * C00;
    result += shc1m1 * (C1 * y);
    result += shc10 * (C1 * z);
    result += shc1p1 * (C1 * x);
    result += shc2m2 * (C2 * x * y);
    result += shc2m1 * (C2 * y * z);
    result += shc20 * (C20 * (3.0 * z * z - 1.0));
    result += shc2p1 * (C2 * x * z);
    result += shc2p2 * (C22 * (x * x - y * y));
    return result;
}

[shader("vertex")]
VertexOutput VsMain(VertexInput i)
{
    VertexOutput o;
    o.position = mul(pushConstant.cameraMatrix, float4(i.position, 1.0));
    o.normal = normalize(i.normal);
    return o;
}

[shader("fragment")]
float4 FsMain(VertexOutput i) : SV_Target
{
    float3 normal = normalize(i.normal);
    normal = float3(-1, -1, 1) * normal.xzy;
    const float3 irradiance = EvaluateIrradianceSH(normal);
    return float4(irradiance, 1.0);
}

SH 係数から irradiance を復元するコードは EvaluateIrradianceSH で計算しています。ここでは単純に SH 係数と対応する SH の値を計算しているだけです。

引数としてここでは normal を渡しています。今回は z-up 座標系の SH で積分したので、座標系を合わせるためちょっと変換しています。これ自体はどの座標系が使用されているかなどによって変わるのであまり気にしなくて大丈夫です。

    normal = float3(-1, -1, 1) * normal.xzy;
    const float3 irradiance = EvaluateIrradianceSH(normal);

これによる結果は以下に示します。Skybox に使用する HDRI は Poly Haven のものを使用しています。

Aarfontein Dusk

Sunset Meadow

Ground truth との比較

どの程度 SH 近似が本来の irradiance と近いかモンテカルロ法で計算した environment irradiance map の ground truth と比較してみます。以下で示す画像は縦を緯度、横を経度として irradiance $E(n)$ をマッピングしたものです。

SH 近似との結果と比較してみると、だいたいの雰囲気が近いことが分かります。SH による近似はかなり精度が良いことが読み取れます。

Ground truth

SH 近似の結果

Ground truth

SH 近似の結果

リンギング

太陽光のような指向性が強い HDRI で実行してみると反対側の色にリング状のアーティファクトが出ることが知られています。例えば次のような HDRI では反対側に不自然に明るい部分が現れます。

これはリンギング（ringing）と呼ばれます。リンギングは指向性に対する SH の精度が足りていないことによるものです。場合によっては負の値が生じることもあります。

リンギングの対策はいくつかあります。ここでは深くは触れませんが、シンプルなものとしては単純に入射光をクランプして指向性を失わせることがあります。もちろんこれは正しい結果を導きませんが、太陽光と HDRI を分離したい場合などには使うことができると思います。（用途次第だと思います）

終わりに

記事を見て頂きありがとうございました。前回から結構時間が経ってしまいましたが、ようやく続きを出すことができました。

実際にやってみると実装はシンプルではありますが、ゲームエンジンでよく見るライティングが自作の viewer 上で実現できるので楽しかったです。

この記事が理解の助けになれば幸いです。

Appendix

実SH

$\ell = 0$

$$
Y_{0,0} = \frac{1}{2} \sqrt{\frac{1}{\pi}}
$$

$\ell = 1$

$$
\begin{align}
Y_{1,-1}(x,y,z) &= \sqrt{\frac{3}{4\pi}} y \\
Y_{1,0}(x,y,z) &= \sqrt{\frac{3}{4\pi}} z \\
Y_{1,1}(x,y,z) &= \sqrt{\frac{3}{4\pi}} x \\
\end{align}
$$

$\ell = 2$

$$
\begin{align}
Y_{2,-2}(x,y,z) &=\frac{1}{2}\sqrt{\frac{15}{\pi}} xy \\
Y_{2,-1}(x,y,z) &=\frac{1}{2}\sqrt{\frac{15}{\pi}} yz \\
Y_{2,0}(x,y,z) &=\frac{1}{4}\sqrt{\frac{5}{\pi}} (3z^2 -1) \\
Y_{2,1}(x,y,z) &=\frac{1}{2}\sqrt{\frac{15}{\pi}} xz \\
Y_{2,2}(x,y,z) &=\frac{1}{4}\sqrt{\frac{15}{\pi}} (x^2 – y^2) \\
\end{align}
$$

$\hat{A}_\ell$の解析解

$$
\begin{aligned}
&l = 1 \quad && \hat{A}_1 = \frac{2\pi}{3} \\
&l > 1,\ \text{odd} \quad && \hat{A}_l = 0 \\
&l\ \text{even} \quad && \hat{A}_l = 2\pi \frac{(-1)^{\frac{l}{2}-1}}{(l+2)(l-1)}
\left[ \frac{l!}{2^l \left( \left( \frac{l}{2} \right)! \right)^2} \right]
\end{aligned}
$$