GPU-Based Importance Sampling

要約

IBLにてレンダリング方程式の積分を近似するために準モンテカルロ法を用いる
現実的な反復回数で真の値に近づけるために、鏡面反射率に応じたImportance Samplingを行う

レンダリング方程式

レンダリングの目的は、ピクセルごとに、ある方向 $u$ に対して、照射された光 $L_i(u)$ がカメラ方向 $v$ に向かって｜双方向反射率分布関数《Bidirectional Reflectance Distribution Function》（BRDF） $f(u,v)$ を持つ面で反射したときの光を、半球 $H$ 上で積分することに相当する。

L_o(v) = \int_H L_i(u) f(u, v) (n \cdot u) du

Image-Based Lighting（IBL）では、この入射光 $L_i(u)$ を環境マップで近似するが、真面目に数値積分法でやろうとすると、テクスチャサンプリングが何百回と必要になり、リアルタイムレンダリングとして重くなりすぎる。

モンテカルロ求積法

｜モンテカルロ求積法《Monte Carlo Quadrature》はランダムな入力値から計算される出力値の相加平均で積分を近似する手法である。理論上、無限回繰り返すと真の値と同値になるが、現実的な有限回数では誤差が大きい。そこで、積分を上手く近似するような最も重要なサンプルをいくつか選択することで誤差を減らすことを考える。

Importance Sampling

光沢のある材質の場合、反射光の多くが鏡面反射方向に反射するため、その周辺に重点を置くのは理に適う。そこで、確率密度関数（PDF）を定義し、鏡面反射率に応じた重み付けを行うことを考える。具体的には以下を用いる。

L_o(u) \approx \frac{1}{N} \sum_{k=1}^N \frac{L_i(u_k) f(u_k, v) \cos \theta_{u_k}}{p(u_k, v)}

ここで、 $p(u_k, v)$ は光源方向 $u_k$ から入射した光が視線方向 $v$ へ反射する確率を表す。

BRDFと入射光の積を用いたサンプリングは正確であるが極めて複雑である。一方、BRDFのみを用いたサンプリングは比較的簡単であり、フィルタリングと一緒に使うことで十分に良い画が得られる。

材質関数のサンプリング

累積分布関数（CDF）の逆関数に一様乱数を入力すると、その分布に従うサンプルが出力として得られることを利用する¹。つまり、サンプルの方向 $(\theta_s, \phi_s)$ を入力とするBRDFを正規化してPDFとし、そのCDFの逆関数を求め、そこに一様乱数のペア $(\xi_1, \xi_2)$ を入力することでサンプル方向 $(\theta_s, \phi_s)$ を求める。

例えば、Phong BRDFのPDFは以下のようになる。

p(\theta_s, \phi_s) = \frac{\cos^n \theta_s \sin \theta_s}{\int_0^{2\pi} \int_0^{\pi/2} \cos^n \theta \sin \theta d\theta d\phi} = \frac{(n+1)}{2\pi} \cos^n \theta_s \sin \theta_s

また、 $p(\theta_s, \phi_s) = p(\phi_s | \theta_s) p(\theta_s)$ から以下が求まる。

\begin{align} p(\theta_s) &= \int_0^{2\pi} p(\theta_s, \phi_s) d\phi_s = (n + 1) \cos^n \theta_s \sin \theta_s \\ p(\phi_s | \theta_s) &= \frac{p(\theta_s, \phi_s)}{p(\theta_s)} = \frac{1}{2\pi} \end{align}

従って、これらのCDFの逆関数から以下が求まる。

\theta_s = \cos^{-1} \left( \xi_1^{\frac{1}{n+1}} \right), \phi_s = 2 \pi \xi_2

ここで、 $n$ は表面のshininessを表す。

Quasirandom Low-Discrepancy Sequences

擬似乱数を用いたサンプリングはうまく分散せずに近似の正確性が損なわれるので、分布の均質性が本質的に保証されているdiscrepancy²の低い｜準乱数列《quasirandom sequences》 で代用する。ここでは、 $N$ 個の二次元のHammersley列を用いる。

x_i = \left( \frac{i}{N}, \Phi_2(i) \right)

ここで、 $\Phi_2(i)$ はradical inverseの関数を表す。これは、2進数表現で小数点を境に $i$ を反転させたような値で、以下のように表される。

$i$ (10進数)	$i$ (2進数)	$\Phi_2(i)$ (2進数)	$\Phi_2(i)$ (10進数)
$1$	$01_2$	$0.10_2$	$0.5$
$2$	$10_2$	$0.01_2$	$0.25$
$3$	$11_2$	$0.11_2$	$0.65$
$\vdots$	$\vdots$	$\vdots$	$\vdots$

MIPMAPを用いたサンプリング

サンプル数が少ないことで起こるエイリアシングを減らすため、MIPMAPを積分の近似として用いる。つまり、PDFの値が小さいほど、サンプル同士が異なる方向を向くようになるため、ピクセルあたりに含まれるサンプルのより多い低詳細度のMIPMAPレベルが選択されるようにする。そのMIPMAPレベルは以下のように計算される。

l = \max \left( \frac{1}{2} \log_2 \frac{\Omega_s}{\Omega_p}, 0 \right) = \max \left( \frac{1}{2} \log_2 \left( \frac{w \cdot h}{N} \right) - \frac{1}{2} \log_2 (p(u,v) d(u)), 0 \right)

ここで、 $\Omega_s = \frac{1}{N \cdot p(u, v)}$ はサンプルに対する立体角を表し、 $\Omega_p = \frac{d(u)}{w \cdot h}$ はレベル0のピクセルに対する立体角を表す。また、 $d(u)$ はマッピングの歪みを表し、環境マップの形式によって異なる。

マッピングと歪み

キューブマップは各面ごとにMIPMAPを生成するため、つなぎ目でアーティファクトが出やすい。そこで、ハードウェアのMIPMAP生成を活用しつつ、つなぎ目を少なくするため、環境マップにdual-paraboloidを用いる。そして、さらにつなぎ目のアーティファクトを減らすため、つなぎ目付近が重複して格納されるようなスケーリングファクター $b$ を導入する（ここでは $b = 1.2$ とする）。このとき、マッピングの歪みは $d(u) = 4b^2 (|z| + 1)^2$ となる。

擬似コード

float4 GPUImportanceSampling(
    float3 viewing : TEXCOORD1
    uniform sampler2D env
) : COLOR {
    float4 c = 0;
    // Sample loop
    for (int k=0; k < N; k++) {
        float2 xi = hammersley_seq(k);
        float3 u = sample_material(xi);
        float pdf = p(u, viewing);
        float lod = compute_lod(u, pdf);
        float2 texcoord = convert_to_dual_paraboloid(u);
        float3 L = tex2Dlod(env, float4(texcoord, lod.xy));
        c += L * f(u, viewing) / pdf;
    }
    return c / N;
}

参考文献

Footnotes

訳注：これを逆関数法（Inverse transform sampling）という。 ↩
訳注：一様分布がどれだけ理想的かを示す指標のようなもの。discrepancyが低いとサンプルが空間的に程よくバラけている。 ↩