こんにちわ、Pocolです。
今日は…
[Xie 2021] Tiantian Xie, “Cache and Bandwidth Aware Real-Time Subsurface Scattering”, University of Maryland, 2021.
を読んでみようと思います。
いつもながら誤字・誤訳があるかと思いますので，ご指摘頂ける場合は正しい翻訳と共に指摘していただけると幸いです。
Chapter 1: Introduction
光が表面に当たると、サーフェイスで跳ね返り、別の場所で跳ね返ることで、柔らかな見た目が生まれます。この効果を表面下散乱といいます。現実の世界では、ろうそく、果物、翡翠などの物体で、驚くほど楽しい表面下散乱の効果があります（図1.2）。しかし、リアルタイムレンダリング、特にゲームにおいて最も重要なのは、フォトリアルな人肌のレンダリング（図1.1および図1.2に示す）であり、人工的なキャラクターが人間に近いが近すぎると恐ろしいことになる不気味の谷に陥らないことです。しかし、フォトリアリスティックレンダリングでは、現実世界での光子の蓄積をシミュレートするために、ピクセルごとに多数のランダムな光子を平均化するモンテカルロサンプリングという高価な数値計算手法が必要です。ゲームの予告編だけでなく、高フレームレートのゲームプレイにおいても、そのようなキャラクターを高い臨場感でレンダリングする技術が切望されています。この制約により、どんなに複雑な散乱であっても、表面下散乱パス全体にはミリ秒、あるいはサブミリ秒の時間しか残されていません。時折、1フレームでもパフォーマンスが低下すると、没入感が損なわれる可能性があります。
表面下散乱をフォトリアリスティックに表現するには、確率的モンテカルロサンプリングに大きく依存することになります。物理ベースのレンダリングをリアルタイムで行うことは、キャッシュアーキテクチャの設計上、インコヒーレントキャッシュアクセスが多く、フレーム内のモンテカルロサンプルの数により高い帯域幅を必要とするため、さらに困難な問題です。
現代および次世代のハードウェアで、高品質なフォトリアリスティック・リアルタイムレンダリングを実現するために、どのような新しい技術を用いることができるでしょうか。これが、本論文の主要な研究課題であり、新たな解決策を提供するために探求します。
1.1 Thesis Statement
モンテカルロサンプリングに基づくフォトリアルな表面下散乱は、現代のGPUキャッシュアーキテクチャにおける非干渉性のキャッシュアクセスとフレーム内の高い帯域幅の要求により、リアルタイムレンダリングでは高コストとなります。数学的に妥当な適応型サンプリングアクセラレーション技術により、サンプル数とメモリ需要を最小化することができます。さらに、スケーラビリティを高めるために、アダプティブサンプリングと分離可能な近似キャッシュのハイブリッドな組み合わせにより、高いフレームレートでの表面下散乱を高品質で実現します。
1.2 Cache and Bandwidth Aware
モンテカルロサンプリングでは、各サンプルアクセスでレンダリングに必要な複雑なマテリアルを評価する場合があります。また、現在のGPUアーキテクチャが持つ時間的・空間的なコヒーレンスが良好であるという前提に沿ったアクセスとは限りません。例えば、スクリーン空間においてより物理的に正しい表面下散乱を持つためには、現在のアーキテクチャ設計に反する放射照度テクスチャを確率的にサンプリングする必要があるだけでなく、ピクセルごとに異なるサブサーフェイスマテリアルをどのようにブレンドすべきかを解決するために、同じ確率的方法で有用テクスチャ（例えば、ピクセル表面下プロファイルおよび/または法線）にもアクセスする必要があります。これは、キャッシュの非干渉性と帯域幅の要求をさらに誇張します。これは、サブサーフェススキャッタリングに限ったことではありません。レイトレーシングを可能にする Turing アーキテクチャ [NVIDIA, 2018] は，キャッシュ・インコヒーレントなシーンコンテキストアクセスの代償として，より物理的に正しいレンダリングを可能にします。交差のハードウェアアクセラレーションはジオメトリアクセスと交差テスト中の計算時間を削減しますが，光と物質の相互作用をシミュレートするために使用されるマテリアルの複雑さは，確率的サンプリングによるリアルタイムレンダリングを妨げます。これは、インコヒーレントなキャッシュの要求と、限られた帯域幅の結果です。本論文では，性能を向上させるために，アルゴリズムレベルでキャッシュと帯域幅を意識した方法を提案しようとするものです。
1.3 Heterogeneous Real-time Rendering
リアルタイムレンダリングを高速化するためのレンダリング技術は数多く存在します。しかし、アルゴリズムレベルでキャッシュとバンド幅の要求を最小化するためのアルゴリズム設計の助けとなるような分類法は存在しません。そこで、リアルタイムレンダリング要求の同質性・異質性に基づき、i) 計算要求, ii) サンプル要求, iii) メモリ要求の3つのカテゴリーに分類しました。この分類法を用いると、リアルタイムレンダリングを高速化する一般的なアルゴリズムの要求とコストを数式で表現することができます。簡単にまとめると、以下のようになります。
1. 計算要求。 演算ユニットコストは、ユニット数、ユニットの種類（例：ラスタライザ、レイトレーシング・ユニット）、アクセラレーション技術の種類（例：Shader levels of detail）、空間再利用要因（可変/適応レートシェーディングによる）に分解することができます。
2. サンプル要求。 サンプルは品質を高めるために使用します。ピクセル間で均質な方法（重点サンプリング、時間的蓄積など）とピクセル間で不均質な方法（適応的サンプリングなど）の両方を用いて、サンプル需要を低減することができます。
3. メモリ要求。 メモリアクセスコストは、プロシージャル技術、圧縮、ミップマップ、仮想テクスチャ、ジオメトリレベルオブディテールによって削減することができます。
分類法の中で、表面下散乱をリアルタイムで効率的に行うためのアルゴリズムを設計しています。この分類法は、キャッシュとバンド幅の観点から設計されていることに注意してください。詳細は第3章をご参照ください。
1.4 Adaptive Sampling
アダプティブサンプリングは、レンダリング効率を上げるための効率的な手法であり、目標とする分散があれば、より少ないレンダリング時間（つまり、サンプル数）で済むようになります。このアルゴリズムは、サンプル数と分散の関係についての統計的定式化に基づいています（つまり、サンプル数が2倍になると分散が半分になる）。その効率性は、分散を調査するために、いくつかのパイロットサンプルを早期に終了させることによってもたらされます。しかし、パイロットサンプルの数が少なすぎると、分散は不正確に推定されます。パイロットサンプル数が大きすぎると、効率が悪くなります。我々は、サンプル履歴を用いた近似関数を使用するリアルタイムレンダリングのためにこの問題に取り組みますが、それでもこの履歴と対応する計算量は\(O(1)\)です。これにより、リアルタイムレンダリングにおけるアダプティブサンプリングの効率が向上します。
このアルゴリズムでは、履歴を利用して効率を向上させながらも、大域的時間蓄積アルゴリズムの実際の実装に制約を与えません。標準的なTAA（Temporal Anti-Aliasing）アルゴリズムのバリエーションで動作することを実証します。我々は、コンピュータグラフィックスにおけるこのグローバルな蓄積のための新しい深層学習アルゴリズムを寄稿していませんが、我々のアルゴリズムは、グローバルな時間的蓄積アルゴリズムとして、スーパーサンプリングのために設計された既存の一般的な深層学習技術、深層学習スーパーサンプリング（DLSS）でも動作することを実証しています。
空間分散を推定するためにコンテキストの時間的調査に基づいて分散推定を行っているため、時間的分散が過大推定につながる可能性があります。本論文では、時間的に安定な適応的サンプリングを実現するために、時間的分散をできるだけ除去するための制御変数を適用しています。オンライン推定を実現するために、簡単でありながら数学的根拠のある新しい手法である指数関数的移動共分散行列をオンライン共分散推定に用いることを提案します。
1.5 Subsurface Scattering
物理ベースのレンダリングによる表面下散乱のポテンシャルを最大限に発揮するために、バーリーの正規化拡散反射プロファイル[Christensen and Burley, 2015]に基づく我々の提案する適応型サンプリング技術でモンテカルロサンプリングを追求します。効率的にするために、重要度サンプリング、および層別サンプリングと低不一致シーケンスの両方を含む、さらにいくつかの技術が利用され拡張されています。
重点サンプリングでは、まず半径重点サンプリング関数を簡略化し、効率的にサンプルを生成します。次に、フィルタ付き重点サンプリング[Křivánek and Colbert, 2008]を適用し、サンプルあたりのサンプル情報を増加させ、帯域幅の要求を減らす。このフィルタ付き重点サンプリングアルゴリズムは、我々の適応サンプリングアルゴリズムと融合し、適応フィルタ付き重点サンプリング(AFIS)として発表されています。これは、表面下散乱は不均質なサンプル要求を持つため、アルゴリズムをサブリニア時間で実行させることで、レンダリング時間をさらに短縮させます。
表面下散乱は、1ピクセル内の直接散乱（拡散）と1ピクセル以上の遠方散乱に分解できるので、それぞれの領域を分離し、層別サンプリングで最終的な値を結合します。分解には2次元サンプリングシーケンス領域での層別サンプリングを適用しています。直接散乱はフラットなピクセル値にアクセスすると仮定しているので、アダプティブサンプリングでは遠方散乱の分散のみを監視しています。
また、アダプティブサンプリングは、すべてのシーンが複雑でない場合にのみ有効であるため、パスである表面下散乱に複数の高速化技術を用いることで、計算負荷を最小化することも追求しています。具体的には、分離可能な表面下散乱とAFISに基づく我々の手法をオンザフライでダイナミックに切り替えるためのフレームワークを提案しました。これにより、メインキャラクターは物理ベースのレンダリング品質をフルに発揮し、他のキャラクターはパフォーマンスと品質のバランスを考慮して分離可能なものを使用するという、高解像度ゲームの可能性を提供します。
1.6 Outline
第2章では、モンテカルロ積分とリアルタイムの表面下散乱の背景と関連する研究を紹介します。
第3章では、キャッシュと帯域幅の要求の観点から、非同次のリアルタイム・レンダリング・アルゴリズムの分類学を提供します。
第4章では、モンテカルロサンプリングのためのリアルタイムアダプティブサンプリングアルゴリズムを紹介します。
第5章では、サブサーフェススキャッタリングフレームワークとその他の高度な設計を紹介します。
第6章では、時間的に安定な適応サンプリングを行うための制御変量に基づくオンラインアルゴリズムを紹介します。
第7章では、表面下散乱の実装と結果の詳細について述べます。
第8章では、結論と今後の研究の方向性を示します。
Chapter 2: Related Work and Background
2.1 Monte Carlo Integration
フォトリアリスティックレンダリングにおける関数の積分（表面下散乱、アンビエントオクルージョン、ソフトシャドウ、光沢反射など）は、解析解がないことが多い。\(F = \int f(x) dx\)として積分する与えられた関数\(f(x)\)に対して、モンテカルロ積分法を用いて、次のように解くことができます。
\begin{eqnarray}
F_N = E [f(x)]_N = \frac{1}{N} \sum_{i=1}^N f(X_i) \tag{2.1}
\end{eqnarray}
ここで、確率変数\(X_i\)は一様分布です。\(X_i\) で関数が評価できる限り、積分は解くことができます。しかし、モンテカルロ推定では、削減率\(O(\sqrt{n})\)の誤差をもたらします。つまり、誤差を半分にするためには、4倍のサンプル数が必要なのです。レンダリングにおける主要な研究課題の一つは、最小限のサンプル数で分散を最小化することです。特に、異なるハードウェアプラットフォームでミリ秒レベルの時間予算が与えられ、帯域幅と計算能力が制限されるリアルタイムレンダリングでは、このような課題があります。単一命令複数データ（SIMD）GPUアーキテクチャをフルに活用するためには、性能の低いデバイスではシェーディングユニットあたりのサンプル数がさらに少なく等しくなることにつながります。したがって、許容できるバイアスが発生するとしても、低コストで分散を低減する技術が決定的に必要です。
2.1.1 Variance Reductio
本節では、確率的サンプリングアルゴリズムを高品質なリアルタイムレンダリングに対応させるために、本アルゴリズムが利用するすべての分散低減技術をレビューします。
2.1.1.1 Importance Sampling
分散低減のためのモンテカルロレンダリングで最も著名な方法の一つは、双方向パストレーシングでグラフィックスに最初に導入された重点サンプリング（IS）です [Lafortune and Willems, 1993, Veach and Guibas, 1995]。基本的な考え方は，ある分布に対する関数の期待値は，次のように別の分布から重み付けされた描画によって近似できるというものです [Tokdar and Kass, 2010]。
\begin{eqnarray}
E_p [f(x)] &=& \int f(x) p(x) dx \tag{2.2} \\
&=& \int f(x) \frac{p(x)}{q(x)} q(x) dx \tag{2.3} \\
&=& E_q \left[ \frac{ f(x)p(x) }{ q(x) } \right] \tag{2.4} \\
&=& E_q [ f(x) w(x) ] \tag{2.5}
\end{eqnarray}
ここで、\(w(x)\)、\(p(x)\)は確率密度関数、\(q(x)\)は別の確率密度関数である。\(E_q[\cdots]\)は期待値が\(q(x)\)に関してであることを表し、\(x〜q\)で\(f(x)q(x) \neq 0\)、\(q(x) \gt 0\)です。そして、\(x\)に基づく関数\(f\)と\(w\)の相関サンプリングに対する分散は次のようになります。
\begin{eqnarray}
Var(f \cdot w) &=& E[f^2 w^2] – E[fw]^2 \tag{2.6} \\
&=& Cov(f^2, w^2) + (Var(f) + E[f]^2)(Var(w) + E[w]^2) \tag{2.7} \\
& & – ( Cov(f, w) + E[f]E[w])^2 \tag{2.8}
\end{eqnarray}
この分散は、\(f\)と\(w\)の間に独立な仮定を持てないため、これ以上小さくすることはできません。しかし、\(f(x)=c\)が定数として知られている場合、次のようになります。
\begin{eqnarray}
Var(f \cdot w) &=& 0 + c^2 (Var(w) + E[w]^2) – c^2 E[w]^2 \tag{2.9} \\
&=& c^2 Var(w) \tag{2.10}
\end{eqnarray}
このような状況下では、\(q\)に対する\(w(x)\)の分散が小さいほど、全体の分散も小さくなり、領域\(x\)にわたって\(q(x)\)が\(p(x)\)に近似するほど、重点サンプリング技術は効率的であることを意味します。
レンダリング、特にリアルタイムレンダリングでは、ターゲット密度関数\(p(x)\)が既知、あるいはCDFが既知であることが多い。前者の場合、\(q(x)=p(x)\)とすることで、ISをフル活用することができる。\(q(x)\)が複雑でリアルタイムにサンプリングできない場合は、軽量な近似関数\(q(x)\)をあらかじめフィッティングしておけばよいです。しかし、後者の場合、解析的な逆解が導けない可能性があるため、CDFを逆解するトライアングルカット解[Heitz, 2020]を適用するか、再度近似を行うことができます。
もう一つの解決策は、複数の密度関数を一緒に採用して\(p(x)\)を近似する多重点サンプリング［Grittmann et al.2019］を使うことです。しかし、照明分布関数のように解析的な密度関数が存在しないこともあります。それらの関数は、データから推定され、構造化重点サンプリング[Agarwal et al., 2003]のようなデータの表現によって記述され、その後、確率変数に基づいて表現をサンプリングすることができます。
リアルタイムレンダリングでは、\(q(x)\)を記述する解析関数があったとしても，確率的にアクセスする必要があるため、オフチップメモリからの頻繁な読み込みが避けられないキャッシュ・インコヒーレンス・クライシスとなり、パフォーマンスを発揮することは容易ではありません。
2.1.1.2 Low-discrepancy Sequence
MC積分のためのサンプルを効率的に生成するためには、サンプリング空間を一様にサンプリングする必要があります。一般に、以下のようになります。
\begin{eqnarray}
E[f(x)] = \underbrace{ \int_{{\mathbf u} \in {\mathcal D}} \frac{f(Q^{-1}({\mathbf u})) p(Q^{-1}({\mathbf u})) }{ q(Q^{-1}({\mathbf u}))} }_{\rm uniform sampling} d{\mathbf u} \tag{2.11}
\end{eqnarray}
積分ワープ\(x = Q^{-1}({\mathbf u})\)とし、\({\mathbf u}\)を一様空間、\({\mathcal D} = [0, 1)^n\)を一様空間とします。\(Q(x) = \int_x q(t)dt\)は累積密度関数です。PDF \(q(Q^{-1}({\mathbf u}))\) は変数変化ヤコビアン行列式である。
\({\mathbf u}\)を均一にサンプリングするために、リアルタイムレンダリングで知覚される誤差を最小にするために、ネイティブのランダムシーケンスではなく、異なる低不一致シーケンス（例えば、Sobol、Halton、Hammersley、\(n-\)Rooksまたは\(R_2\)）が好まれます。シーケンスの不一致度[Phar et al., 2016]は、次のように測定することができます。
\begin{eqnarray}
D_N(B, P) = \underset{b \in B, B \subseteq {\mathcal D}}{\rm sup} \left| \frac{ \# \{ {\mathbf u}_i \in b\} }{N} – V(b) \right| \tag{2.12}
\end{eqnarray}
ここで、\(N\) 個のサンプル \({\mathbf u}_1, \cdots , {\mathbf u}_N\) のシーケンスは、\(b\), \(B = \{[0, v_1] \times \cdots \times [0, v_i] \times \cdots \times [0, v_n] \}\)で\(v_i \in [0, 1)\) 内のシーケンス発生器に従って描画されます。\(V(b)\)は注目量を表し、\(\#\{\cdot\}\)はサンプル数を表します。
リアルタイムレンダリングパイプラインで、時間的蓄積の機能がある場合、最適な低乖離シーケンスは、フレーム内のモンテカルロ積分のための投影\((n-1)\)次元空間でも低乖離である低ディスクレパンシー\(n\)次元数列となります。我々の知る限り、これはまだ文献上未解決の問題である。本論文では、リアルタイムレンダリングのためのサンプリング配列の改良については検討しません。しかし、サンプル生成には、リアルタイムに適した低ディスクレパンシー数列である \(R_2\) 数列を使用する予定です。近傍ピクセル間の相関を緩和するために、リアルタイムに効率の良いハッシュ関数[Jarzynski and Olano, 2020]を適用し、例えばフレーム番号とピクセル座標に基づいてランダムサンプルを生成する。
2.1.1.3 Stratified Sampling
層化サンプリングも分散を減らすのに有効な手法の一つで、ドメインを層状に細分化します。各層は、異なるサンプリング戦略でサンプルを割り当てることが保証されています。そして，修正MC積分において分散を最小化することができます。
\begin{eqnarray}
E[f(x)] = \sum_{i=1}^M \int_{ {\mathbf u} \in {\mathcal D}_i } \frac{ f(Q^{-1}({\mathbf u})) p(Q^{-1}({\mathbf u})) }{ q(Q^{-1}({\mathbf u})) } d{\mathbf u} \tag{2.13}
\end{eqnarray}
ここで、\({\mathcal D} = {\cup}^M_i {\mathcal D}_i\)と各\({\mathcal D}_i\)は層である。この地層境界の細分化は、コンピュートグラフィックスではグリッドとして一定であることが多いです[Phar et al., 2016]。しかし、\(P_i = \frac{ |{\mathcal D}_i | }{ |{\mathcal D}| }\).がある確率密度関数に基づいて割り当てることができます。\(| \cdot |\)は体積を測定します。その場合、\(N\) サンプル MC 推定量は
\begin{eqnarray}
F_N = \sum_{i=1}^M P_i \left[ \frac{1}{N_i} \sum_{j=1}^{N_i} \frac{ f(Q^{-1}({\mathbf u}_j)) p(Q^{-1}({\mathbf u}_j)) }{ q(Q^{-1})({\mathbf u}_j)) } \right], {\mathbf u}_j \in {\mathcal D}_i \tag{2.14}
\end{eqnarray}
ここで、\(N_i = P_i \cdot N\). 層別化により、重要な特徴が全体としてサンプリングされる割合が高くなるようにします。さらに、層別化は、いくつかの層別積分が、少ないサンプル数（例えば、1サンプル）で正確な結果が得られる低頻度であることが分かっている場合、サンプリング効率を向上させます。これは実際に、5.3.1節で述べた散乱の統一による、より効率的なリアルタイム適応サンプリングの設計につながるものです。
2.1.1.4 Control Variates
制御変量法は、モンテカルロ積分を、あるスケーリング係数を持つ既知の積分と、より簡単に推定できる残差の和に分解する分散低減法です [Ripley, 2009]。金融 [Alexader, 1999]、オペレーションズリサーチ [Nelson, 1990, Hesterberg and Nelson, 1998]、コンピュータネットワーク [Lavenberg et al., 1982] など、多くの分野で研究・応用されています。最適な CV 係数の推定が必要です。
積分する関数\(F = \int_{\mathcal D} f(x) dx\)に対して、与えられた領域Dの積分がすでに\(G\)として知られている制御変量関数\(g(x)\)があると仮定し、その関数は次のように書き換えることができます。
\begin{eqnarray}
F = a \cdot G + \int_{\mathcal D} f(x) – a \cdot g(x) dx \tag{2.15}
\end{eqnarray}
ここで、\(a\)は信号\(f(x) – a \cdot g(x)\)ができるだけ平坦になるように関数\(g(x)\)をスケールするためのパラメータです。サンプリングの際、f(x)とg(x)は、結合サンプリングと同じ乱数でサンプリングされます。このようにすると、信号の平坦性により、xの違いは比較的小さくなります。分散は信号内部の傾きによってのみもたらされます。そして，CV推定量は
\begin{eqnarray}
\langle F \rangle = a \cdot G + \langle F – a \cdot G \rangle \tag{2.16}
\end{eqnarray}
\(\langle \cdot \rangle\)は、簡単のためにMCサンプリング期待値を表します。以下として低分散であることを示す必要があります。
\begin{eqnarray}
Var(\langle F \rangle) &=& Var(a \cdot G + \langle F – a \cdot G \rangle) \tag{2.17} \\
&=& Var( \langle F \rangle) + a^2 Var (\langle G \rangle) – 2 a Cov ( \langle F \rangle, \langle G \rangle), \tag{2.18}
\end{eqnarray}
ここで，\(a = Cov(\langle F \rangle, \langle G \rangle) / Var(\langle F \rangle)\) のとき（式 2.18 の微分を 0 にして a を解けば），分散は次のように最小化される。
\begin{eqnarray}
Var(\langle F \rangle) = (1 – Corr ( \langle F \rangle, \langle G \rangle)^2) Var (\langle F \rangle) \tag{2.19}
\end{eqnarray}
ここで、\(Corr(\cdot)\)は相関係数であす。したがって、選択した制御変量関数\(g(x)\)が\(f(x)\)と相関がある限り、分散低減を達成することができます。\(Corr(\langle F \rangle, \langle G \rangle)\)が1または-1になるほど\(g(x)\)が\(f(x)\)に近似していれば、分散はゼロになります。相関がゼロの場合でも、測定された分散は\(Var(\langle F \rangle)\)と全く同じになります。この特性は，分散が推定プロセスで使われるとき，CVを使う非常によい動機を与えてくれます。
コンピュータグラフィックスでは，オフラインレンダリングにおいて，主要部分を分離するためにあらかじめ定義された一定のCV係数を持つアプリケーション（例えば，残差比トラッキング［Novakら，2014］，多重相関サンプリング［Szecsiら，2004］）が多く検討されています．また、ペナルティ付き最小二乗法 [Fan et al., 2006]、回帰 [Rousselle et al., 2016、Kondapaneni et al., 2019]、深層学習 [Muller et al., 2020] によって最適なCV係数を求めることが検討されています。しかし、フレームあたりの時間予算が数ミリ秒しかなく、それらの高価な計算をする余裕がないリアルタイムレンダリングのためのガイダンスやアプリケーションを提供する既存の研究はありません。さらに、ここでの目標は少し異なり、それらの先行技術が静的なシーンで行うように、MC推定量の分散を減らすのではなく、時間的変化による監視分散を減らすことです。監視された分散は、リアルタイムアダプティブサンプリングのサンプル数推定に使用されます。本論文では、アダプティブサンプリングの文脈で、新しいオンライン共分散行列式に基づく、このようなオンラインCV係数推定アルゴリズムを提案します。従来の利用法では、\(f – a \cdot g\)は定数であることが予想されます。しかし、グラフィックスで使われるこの知見は、最適な結果を導き出すものではありません。代わりに時間情報を考慮すると0になることが予想されます。また、簡単に採用できるように、主要部分分離の精神に基づいた定数軽量近似を提供します。詳しくは6章を参照してください。
2.1.1.5 Temporal Reuse
リアルタイムレンダリングは、高い時空間コヒーレンスを示します[Scherzer et al.2012]。このコヒーレンスは，テンポラルアンチエイリアス [Karis, 2014] (TAA) や，コストを償却してレンダリング品質を向上させるためのシェーディングやジオメトリ情報との累積に利用されてきました。文献で扱われている主な問題は，履歴を現在のフレームとブレンドする指数移動平均の使用と，それらが同じオブジェクトからのものかどうかを見分けるヒューリスティックによるブラー，ゴースト，ラグ，フリッカリング［Xiao et al, 2018, Patney et al, 2016, Iglesias-Guitian et al, 2016］です。
また、時間的な情報は他のシステムと統合し、サンプルコストの償却も行っています。最近の研究では、ゴーストやブラーが分散と深度に基づいて検出される場合、TAAの結果をレイトレーシングに置き換えます[Marrs et al., 2018]。時空間分散誘導型フィルタリング（SVGF）法は、ジオメトリ情報と画像空間ウェーブレットフィルタを含み、1ピクセルあたり1サンプルのリアルタイムパストレーシングを実現します［Schied et al.、2017］。Schiedら[2018]は、時間的勾配に基づいて指数的重みを適応的に変更し、SVGFを使用する際に速い光の変化によるシーンの時間的安定性を改善（テンポラルオーバーブラーリングを軽減）します。
本研究では、時間的に蓄積された分散と有効サンプル数によって、確率的サンプリングプロセスをガイドする。これにより、ピクセルあたりの品質要件と時間予算を満たすために、サンプル数を償却することができます。これは、既存のTAAプロセスを変更したり、その実装を制限したりすることなく、TAAプロセスに対してより安定で低分散な入力を提供する一方で、フレーム間でサンプルを結合するためにTAAに依存します。線形モデルと逆共分散行列が分散情報を提供するために再帰的最小二乗法によって適合される既存の研究[Iglesias-Guitian et al., 2016]とは異なり、我々はフレーム間で渡される1つの追加のテクスチャのみを必要とする軽量のアルゴリズムを使用します。
2.1.2 Adaptive Sampling
適応型サンプリングは、ばらつきの大きい領域に多くのサンプルを割り当てるサンプリング手法である[Phar et al.2016]。一様サンプリングよりも高速な収束を可能にします[Dammertz et al., 2010, Moon et al., 2014]。ただし、パイロットサンプルも平均値の推定に使用すると、結果が偏ってしまうという大きな欠点があります[Kirk and Arvo, 1991]。これはまた，リアルタイムレンダリングにおけるアダプティブサンプリングの採用を妨げるジレンマでもあります – ほとんどのリアルタイムレンダラーは，未レンダリングのパイロットサンプルにサンプル数を浪費したくはないのです。このセクションでは，いくつかの適応サンプリングアルゴリズム，バイアスの原因，および適応サンプリングのいくつかの応用についてレビューします。
特に、リアルタイムレンダリングに適した適応型サンプリングアルゴリズムに関心があります。すなわち、適応サンプリングアルゴリズムのオーバーヘッドを犠牲にして、サンプリング効率を上げるようなアルゴリズムであってはなりません。これらのアルゴリズムの中で、我々は統計的分散法が数学的に健全であり、リアルタイムレンダリングに効率的であると判断しています。そこで、その理論と潜在的な偏りについて完全な紹介をします。最後に、将来の潜在的な研究の方向性を示すために、他の潜在的な測定基準について簡単に紹介します。
2.1.2.1 Theory
統計的に最適化された適応サンプリングが文献で提案されている[Lee et al., 1985]。関数\(f(x)\)を領域\({\mathcal D} = R^n\)上で重み関数\(p(x)\)を用いて積分すると、期待値は
\begin{eqnarray}
E_p[f(x)] = \int_{\mathcal D} f(x)p(x)dx, \int_{\mathcal D} p(x)dx = 1 \tag{2.20}
\end{eqnarray}
となります。\(x_1, x_2, \cdots, x_N\) が \(p(x)\) に従って独立に同分布する標本であるとき、モンテカルロ推定量は
\begin{eqnarray}
F_N = \frac{1}{N} \sum_{i=1}^N f(x_i) \tag{2.21}
\end{eqnarray}
となります。大数の強法則より, \({\rm lim}_{N \rightarrow \infty} F_N = E[f(x)] = E[F_N]\)で，そして
\begin{eqnarray}
Var(F_N) &=& E[(F_N – E[f(x)])^2] \tag{2.22} \\
&=& \frac{ E[f^2(x)] – E[f(x)]^2 }{ N } \tag{2.23} \\
&=& \frac{ Var(f(x)) }{ N } \tag{2.24}
\end{eqnarray}
\(N\)の数が増加するにつれて，モンテカルロの分散は低減します。正規理論のもとで分散を所定の閾値以下にするために\(N\)が十分に大きい\(\beta\)確率は
\begin{eqnarray}
P \left( \frac{Ns^2_N}{Var(f(x))} \lt {\chi}^2_{\beta} (N – 1) \right) = \beta \tag{2.25}
\end{eqnarray}
ここで，\(s^2_N = \frac{1}{N}E[(f(x_i) – F_N)^2]\) はデータの標本分散です。なお、カイ二乗は正規理論では次のような形をしています。
\begin{eqnarray}
{\chi}^2(k -1) = \frac{(k-1) \cdot s^2}{\sigma^2} \tag{2.26}
\end{eqnarray}
そして\(k\)は自由度数です。モンテカルロ推定量\(F_N\)の分散に関して、以下のようになる。
\begin{eqnarray}
P \left( Var(F_N) \lt \frac{s^2_N}{ {\chi}^2_{\beta} (N-1) } \right) = 1 – \beta \tag{2.27}
\end{eqnarray}
すると、分散閾値は次のように求めることができます。
\begin{eqnarray}
{Var}_0 (F_N) = \frac{ \sigma^2_w }{ {\chi}^2_{\beta} (N_{\rm max} – 1) } \tag{2.28}
\end{eqnarray}
ここで，\(\sigma^2_w\) は目標の最大分散，\(N_{\rm max}\) は計算バジェットの関係で許される最大サンプル数です。そして、分散の閾値に基づいて最小サンプル数\(N\)を求めることができます。
カイ二乗の再解釈 式2.26のカイ二乗の正規理論の仮定のもとで、データ\(s\)の標準的な導出で\(k\)個の観測値を与え、正規母集団の分散\(\sigma^2\)に達するために必要なサンプル数という、適応的サンプリングの領域での新しい解釈を与えています。我々の提案するリアルタイム適応サンプリングは、この再解釈から式（4.1）に拡張されたものです。
2.1.2.2 Bias
適応的サンプリングは、初期サンプリングで高周波情報を観測できない場合、パーマネントストップが発生します。この早すぎる停止は、バイアスの原因となります。適応サンプリングによる偏りを統計的に説明するために、適応サンプリングアルゴリズム[Kirk and Arvo, 1991]の一般的な形をアルゴリズム1に示す。このアルゴリズムは適応的サンプルのほとんどのケースを捉えている。このアルゴリズムでは，まず，分布のある確率変数\(X\)からサイズ\(n\)のパイロットサンプル集合が引かれる（2行目）。標本集合\(S_n\)の分散指標\(Variantion(\cdot)\)が分散閾値より小さい場合（3行目）、その標本集団は容易とみなされることになります。近傍の値は低頻度です。サンプリング結果はサンプル集合の平均メトリックスとなります（4行目）。そうでなければ、母集団は、近傍値が高頻度であるハードケースとなります。真の平均を計算するために、さらにオラクル演算が行われます（6行目）。
レンダリングにおけるピクセルは、そのほとんどが有限の領域（例えば、対象画素の周囲の補助バッファ）によって決定されるので、確率変数\(X\)は補助バッファの尺度の重み付き和であると仮定して差し支えないと思われます。同様に、それは\(k\)個の異なる値\(v_1, v_2, \cdots, v_k\)を対応する確率\(p_1, p_2, \cdots, p_k\)で総和したものです。
このとき、グランドトゥルース平均は次のように表すことができます。
\begin{eqnarray}
v = \sum_{i=1}^k p_i v_i \tag{2.29}
\end{eqnarray}
ここで、重みは\(\sum_{i=1}^k p_i = 1\)で正規化されます。アルゴリズム1に関して、低周波、高周波に応じた条件付き期待値に基づく期待値を次のように書き換えることができます。
\begin{eqnarray}
E[\mu] = E[\mu | low] \cdot P[low] + E[\mu | high] \cdot P[high] \tag{2.30}
\end{eqnarray}
ただし，\(E[\mu | high] = v\)で\(P(high) = 1 -P(low)\)です。任意の標本集合\(S_n\)は、\(\sum_{i=1}^k n_i = n\)を含む\(k\)個のタプル\((n_1, n_2, \cdots, n_k)\)で作ることができ、\(n_i\)個の標本は値\(v_i\)であるから、多項分布に基づくk-タプルの確率は次のようになります。
\begin{eqnarray}
P[n_1, n_2, \cdots, n_k] = \frac{n!}{n_1! \cdots n_k!} \prod_{i=1}^{k} p_i \tag{2.31}
\end{eqnarray}
低周波の場合、\(\epsilon\)が十分に小さければ、\(n\)個のサンプルはすべて同じ値になる可能性があります。それは次のようになります。
\begin{eqnarray}
P[low] &=& \sum_{i=1}^k p^n_i \tag{2.32} \\
E[\mu | low] &=& \frac{ \sum_{i=1}^k p^n_i v_i }{ P[low] } \tag{2.33}
\end{eqnarray}
適応サンプリング結果関数は次のように簡略化されます。
\begin{eqnarray}
E[\mu] = v + \sum_{i=1}^k p^n_i (v_i – v) \tag{2.34}
\end{eqnarray}
ここで、右の部分はバイアスです。\(n\)が増加すると、バイアスは一貫して減少する。バイアスを最小にするためには、最初のパイロットサンプルサイズ\(n\)を十分に大きくする必要があります。サブピクセル・アンチエイリアスのようなアプリケーションでは、数個のパイロットサンプルで十分な場合があります。しかし，パーセンテージクローザーソフトシャドウ，アンビエントオクルージョン，またはサブサーフェススキャッタリングのような複雑なアプリケーションでは，さらに大きな\(n\)が必要になり，適応サンプリングは不適切になります。本論文では、リアルタイムアプリケーションのために、この問題に取り組むことを試みます。これは、我々の適応サンプリング処理が、各フレームのシェーディング結果ではなく、サンプルカウントを決定するためです（4章）。
2.1.2.3 Metrics
ここでは、適応的サンプリング分散メトリクスを紹介する。Zwickerら[2015]は、メトリクスの種類を「a priori」と「a posteriori」に分けて、最新の適応型サンプリングと再構成アルゴリズムを体系的にレビューしています。
「推測メトリクス」 推測メトリクスは、光輸送方程式の局所的な分析に依存している。例えば、ローカルジオメトリ、BRDFモデル、深度[Soler et al., 2009]、ローカルライトフィールドによって記述することができます。推定後，それらはしばしばサンプル数を決定するために直接使用することができます。具体的には，1）フーリエ領域で周波数解析を行い［Durand et al., 2005］，その周波数情報を1つのサンプルに合成すると，ナイキストのサンプリング定理とローカルバンドリミットに基づいて，局所サンプリング密度を決定することが可能です。また、同じ考え方はウェーブレット領域にも展開される。Bagherら[2012]は、サンプリングレートを推定するために帯域幅バッファを抽出しますが、実際のサンプリングの前に少なくとも16 sppのコストを要します。このスペクトルは、高次元空間の共分散行列によって記述することもできます[Belcour et al.］。 2) 微分解析（1次、2次微分など）も、勾配が高いほど表面が複雑になり、より多くのサンプルが必要になるという基本的な考えで行うことができます。Ramamoorthiら[2007]は、1次勾配の大きさを利用して、一様格子の階層的サンプリングを行っています。勾配が高いほど、次のレベルではより多くのサンプルが抽出されます。
「帰納メトリクス」 これらの方法は、光輸送方程式の実際のモンテカルロサンプリングに依存し、ピクセルまたは領域についてより多くのサンプルがまだ必要であるかどうかを統計的に決定するものである。共通の特徴として、コントラスト [Mitchell, 1987]、知覚的メトリック [Bolin and Meyer, 1998]、MSEによる分散、相互情報、カイ二乗があります。メトリックを決定するために、空間情報が使用されます。例えば、MSEでは、分散を推定するために、近傍のピクセルにアクセスすることが必要です。Bolin と Meyer [1998]は，ウェーブレット空間における階層的な適応サンプリング手法を提案しました。これは，中間誤差を格納した四分木を走査した後，グローバルに最も知覚誤差の大きい位置にサンプルを配置するものです。
他のアルゴリズムでは、「推測」と「帰納」の両方の指標を利用するものもあります。例えば、非一様な単一散乱媒体に対して、Renら[2008]は勾配に基づく補間を用いて、少ないサンプルセットでボリューム内の入力放射輝度を再構成している。しかし、より高密度なサンプリングが必要とされる鋭い陰影領域が存在します。適応的なサンプリングを行うために、彼らは、基底方向に沿った周囲6点の放射輝度値の数値勾配と、有効球内の近隣のサンプリングされた放射輝度間の実際の放射輝度差の積に基づく誤差行列を提案しました。そして、アルゴリズム1の再帰的バージョン（Oracle関数は関数そのものに置き換えられ、新しいサンプルが\(S_N\)に連結される）が、所定の再帰回数に達するまで実行されます。
2.1.2.4 Adaptive Rendering
適応型サンプリングの最近の応用例としては、適応型レンダリングが挙げられる。サンプリング密度やフィルタリングを調整することで、レンダリング効率を向上させることができます。疎な線形モデルと最適なサンプリングウィンドウは，隣接するピクセルを最適に再構成するために再帰的に更新されます[Moon et al. 2015］。誤差推定器によって，高分散領域へのサンプリングが導かれます。適応的サンプリングアルゴリズムは，最適化関数を解くために，例えばパストレーシングのための勾配領域[Manzi et al.，2016]のような大きなメモリと計算能力を必要とすることが多いです。一般的に、2つのパスが必要とされます[Moon et al., 2014]。最初のフェーズでは、誤差を収集するために少数の均一なサンプルが割り当てられます。2回目のパスでは，1回目のパスでのエラーメトリクスに基づき，追加のサンプルが割り当てられます。適応サンプリングの拡張版として，ある変動基準を満たすまでサンプリングを継続する反復適応サンプリングがあります。1パスで済むような工夫を求めます。フレーム間で1枚のテクスチャを追加保存するだけでよく、計算オーバーヘッドも少ないシングルパス適応型サンプリングモデルを提案します。提案する1パスモデルは、動的な光や動きの変化を拾って、効果的にサンプルを調整することができます。最近，NVIDIAの研究者は，特徴バッファ（例えば，法線，深度，ファーストヒット時のアルベド）と前フレームの投影レンダリングに適応的サンプリングメトリクスとしてニューラルネットワーク（U-Net）を適用し，同様にシングルパスでサンプルカウントマップを推定しました［Hasselgren et al, 2020, Hofmann et al, 2021］。しかし，あまりにも高価であり，この適応的なサンプリングタスクには適していません。
2.2 Subsurface Scattering
このセクションでは、まず、リアルタイムレンダリングのための物理ベースのサブサーフェス散乱に関する研究文献をレビューします。次に，様々なモデルやリアルタイムアクセラレーション技術に関する背景を紹介します。
表面下散乱のための初期の物理ベースのレンダリングアルゴリズムは、放射伝達方程式を解くためにモンテカルロパストレーシングを使用していました。BSSRDF（Bidirectional Scattering-Surface Reflection Distribution Function）［Nicodemusら，1977］を導入してサーフェイス領域へのサブサーフェス積分を単純化しても，これらのアルゴリズムは当時のハードウェアでは画像生成に数時間かかっていました．Jensen は、BSSRDF の光拡散モデルとして、ダイポール近似 [Jensen et al., 2001] とマルチダイポール近似 [Donner and Jensen, 2005] を導入しました。より単純なダイポール法は、高散乱媒質における等方的な表面下散乱のみを考慮しますが、拡散反射率プロファイルを導入し、コンパクトな表現と効率的な評価を可能にしました。
拡散プロファイルを使用する場合、リアルタイムレンダリングの複雑さを軽減するために、オンラインサンプリングと事前積分という2つの異なる方向性があります。オンラインサンプリングは，表面下散乱がアーティストフレンドリーカーネル[Borshukov and Lewis, 2005]またはダイポールベースのガウスカーネルの加重和によって近似できることを前提に設計されています [d’Eon et al., 2007]。サンプル数は，ディファードレンダリング後にサブサーフェイスオブジェクトごとにテクスチャ空間［d’Eon and Luebke, 2007］からスクリーン空間［Jimenez et al, 2009, 2015］にシフトすることでさらに削減されます。対照的に，事前積分法は，光勾配変化があるときにサブサーフェス散乱が見えるという仮定[Penner and Borshukov, 2011]と影への距離が正しく再構築できることに基づいて，ダイポール近似に基づくオンラインサンプリングをテクスチャルックアップにオフロードします。
物理ベースのサブサーフェイスレンダリングにおける最近の進歩は、すでに近似であるダイポールモデルを当てはめる代わりに、モンテカルロ結果に近似を当てはめます[Christensen and Burley, 2015]。この手法は、最近、事前に計算されたサンプリングでUnityゲームエンジンに組み込まれました[Golubev, 2018]。
2.2.1 Real-time Model
サーフェイス点\(p\)において、方向\(\omega_o\)に対する出射放射輝度は次のようになります。
\begin{eqnarray}
L_o(p, \omega_o) = \int_{\partial {\mathbf \Omega}} \int_{{\mathbf S}^2} L_i(q, \omega_i) S(q, \omega_i, p, \omega_o) d \omega_i dq, (q \in {\partial} {\mathbf \Omega} \tag{2.35}
\end{eqnarray}
ここで、\(L_i\)は点\(q\)における方向\(\omega_i\)からの入射放射輝度です。\(\omega_i\)は半球\({\mathbf S}^2\)のすべての入射方向について積分し、\(q\)は局所サーフェイスパッチ\(\partial {\mathbf \Omega}\)について積分します。\(S\)は双方向散乱表面反射率分布関数（BSSRDF）である。\(S\)は、方向\(\omega_i\)からの入射光束が\(q\)にあるとき、方向\(\omega_o\)の\(p\)における反射光束の比率を表しています。
\begin{eqnarray}
S(q, \omega_i, p, \omega_o) = \frac{dL_o(p, \omega_o)}{d\Phi_i (q, \omega_i)} [m^{-2} \cdot sr^{-1}] \tag{2.36}
\end{eqnarray}
式2.36は、均質な半無限平面媒質中の半径方向に対称な \(S\) を仮定して、次のように簡略化されることが多いです。
\begin{eqnarray}
S(q, \omega_i, p, \omega_o) = C F_t(q, \omega_i) R(r_q) F_t(p, \omega_o) \tag{2.37}
\end{eqnarray}
ここで、\(r_q = || p – q ||\)で，\(F_t(\cdot)\)は方向性フレネル透過項、\(R(\cdot)\)は拡散反射率プロファイル, \(C\)は定数項です。そして、\(p\)における任意の放射束の反射方向\(\omega\)に対して、サーフェイス点\(p\)において対称な拡散プロファイルを持つ簡略化した表面下散乱関数は、以下のようになります。
\begin{eqnarray}
L_o(p, \omega) &=& \int_{\partial {\mathbf \Omega}} R(r_q) \cdot b (p, q, \omega) \cdot dq \tag{2.38} \\
&=& b (p, q, \omega) = \int_{{\mathbf S}^2} C F_t(p, \omega) F_t(q, \omega_i) L_i(q, \omega_i) \langle \omega_i, n_q \rangle d\omega_i \tag{2.39}
\end{eqnarray}
この式を使用するには、例えば、遅延レンダリングで、内側積分は、\(b\)としてすべてのライトの放射照度テクスチャに事前に計算することができます。リアルタイムレンダリングの文献では、拡散反射プロファイルがどのようにモデル化されるかについて、2つの手法があります。
2.2.1.1 Dipole Profile
Grosjean [1959]以降のほとんどの研究 [Jensen et al., 2001, Donner and Jensen, 2005, d’Eon and Irving, 2011, Hable et al., 2013] は、光（または中性子）を、強度の減少、単一散乱、多重散乱項に分解し、次のように記述しています。
\begin{eqnarray}
S = S^{(0)} + S^{(1)} + S_d \tag{2.40}
\end{eqnarray}
この式はノイマン展開に基づき、別々に処理します。削減された強度透過率\(S^{(0)}\)は散乱のない直接放射輝度であり、Lambert-Beerの法則に基づく標準的な透過率アルゴリズムで実装することが可能です。\(S^{(1)}\)は、メディア内のすべての単一散乱イベントの積分です[Jensen et al.、2001]。多重散乱項\(S_d\)は、残りのすべての散乱事象をモデル化し、ダイポール拡散近似で簡略化されています。古典的なダイポールに基づいて対応する拡散反射率プロファイル\(R_d\)は次のとおりです。
\begin{eqnarray}
R_d(r) = \frac{ \alpha’ z_r(1 + \sigma_{tr} d_r) e^{-\sigma_{tr} d_r} }{ 4 \pi d^3_r } – \frac{ \alpha’ z_v (1 L \sigma_{tr} dv) e^{-\sigma_{tr}d_v} }{ 4 \pi d^3_v } \tag{2.41}
\end{eqnarray}
ここで、\(\alpha’ = \frac{{\sigma_s}’}{{\sigma_t}’}\)は減弱されたアルベドです。\({\sigma_s}’\)と\({\sigma_t}’\)は減弱された散乱係数と減弱された消滅係数で\({\sigma_t}’ = \sigma_a + {\sigma_s}’\)です。\(\sigma_a\)は吸収係数です。\(\sigma_{tr} = \sqrt{3 \sigma_a {\sigma_t}’}\)は有効輸送係数です。\(z_r = \frac{1}{{\sigma_t}’}\)で，\(z_v = \frac{ 1 + 4(1 + F_{dr}) / (3(1-F_{dr})) }{ \sigma + t’ }\) は\(z=0\)を基準とした正負のソース位置であり，\(d_r\)，\(d_v\)は点\(q\)からの対応するソースの距離を\(d = \sqrt{r^2+z^2}\)としたものです。\(F_{dr}\)は拡散フレネル反射率で、多項式展開[Egan et al., 1973]により次のように近似されます。
\begin{eqnarray}
F_{dr} \approx \begin{cases}
0.4399 + \frac{0.7099}{\eta} – \frac{0.3319}{\eta^2} + \frac{0.0636}{\eta^3}, & \eta \lt 1 \\
-\frac{1.4399}{\eta^2} + \frac{0.7099}{\eta} + 0.6681 + \frac{0.0636}{\eta^3}, & \eta \gt 1
\end{cases} \tag{2.42}
\end{eqnarray}
ここで，\(\eta\)は屈折率の比率です。
コンピュータグラフィックスの文献には、より正確なダイポールモデル（マルチポールなど）がありますが、リアルタイムレンダリングの範囲では深い内容をカバーせず、Burleyのプロファイルに主眼をおいています。
2.2.1.2 Burley’s Profile
散乱事象を個別に近似する代わりに，ChristensenとBurley [2015]は，すべての散乱項を含む経験的MCシミュレーションデータに基づいて拡散反射率プロファイル\(R(\cdot)\)を直接近似しています。\(R(\cdot)\)は、次のように距離\(r\)に関して2つの指数関数の和でうまく近似することができます。
\begin{eqnarray}
R(r) = A \frac{ e^{-r/d} + e^{-r/(3d)} }{ 8\pi dr } \tag{2.43}
\end{eqnarray}
ここで、\(A\) はサーフェイスアルベドで、\(A=\int_0^{\infty} R(r) 2 \pi r dr\) [Jensen and Buhler, 2005] とすると、Burley の正規化プロファイルとすることができます。\(d\)項は設定\(\Theta\)に基づく自由行程長パラメータ\(\ell\)に適合します。
\begin{eqnarray}
d =
\begin{cases}
\ell / (1.85 – A + 7 |A-0.8|^3) & \Theta = \Theta_1 \\
\ell / (1.9 – A + 3.5(A – 0.8)^2) & \Theta = \Theta_2 \\
\ell / (3.5 + 100(A – 0.33)^4) & \Theta = \Theta_3
\end{cases} \tag{2.44}
\end{eqnarray}
ここで、以下のようになります：
- サーチライト設定(\(\Theta = \Theta_1\)). 光はサーフェイスに対して垂直に体積に入射する。\(\ell\)は体積平均自由経路です。
- ディフューズサーフェイストランスミッション(\(\Theta = \Theta_2\)). \(\ell\)は理想的な拡散透過後のラフマテリアルの体積平均自由経路です。
- パラメータとしてのディフューズ平均自由経路(\(\Theta = \Theta_3\)). \(\ell\)はサーフェイス上の拡散平均自由経路であり、\(\ell_d\)と表記されます。
一般的なゲーム用途では、よりアーティストフレンドリーな\(\Theta_3\)が選ばれています。この選択により、現実世界[Jensen et al., 2001]からBurleyのモデルにフィットしたマテリアルパラメータを以下の式で変換することができます[Christensen and Burley, 2015]。
\begin{eqnarray}
\ell_d &=& \frac{1}{\sigma_{tr}} \tag{2.45} \\
\sigma_{tr} &=& \sqrt{ \frac{\sigma_a}{D} } \tag{2.46} \\
D &=& \frac{\sigma_t + \sigma_a}{3 {\sigma_t}^2} \tag{2.47}
\end{eqnarray}
ここで、\(D\) は拡散係数です。\(\sigma_t = \sigma_a + \sigma_s\)は消衰係数、\(\sigma_s\)は散乱係数です。\(\sigma_s’ = \sigma_s(1 – g)\)は散乱角の平均余弦である\(g\)の影響を受ける縮小散乱係数で、光が前方\((g \lt 0)\)または後方\(g \gt 0)\)に散乱する様子をパラメータ化したものです。\(g=0\)の場合、光は各散乱イベントで一様に散乱されます。一般的なマテリアルの拡散平均自由経路の計算結果[Jensen et al.，2001]は、\(g=0\)の場合、表C.1に示されています。
2.2.2 Real-time Acceleration Techniques
GPU で効率的に動作させ、よりフォトリアルにするために、表面下散乱はテクスチャ空間へ、そしてスクリーン空間へと移行してきました。リアルタイムレンダリングエンジンでは、ガウスブラー、前積分、セパラブルなどの高速化技術が多く採用されています。本論文では、Adaptive Filtered importance sampling という新しい適応的サンプリング手法に基づき、より高いレンダリング品質でサブサーフェスを高速化する手法を提案します。また，この技術は他のレンダリングタスクにも適用可能です。
2.2.2.1 Gaussian Blur
リアルタイムのリアルな表面下散乱は，Borshukov と Lewis [2005] によって，映画「マトリックス・リローデッド」の人間の顔のレンダリングのために初めて研究されました。彼らは、パストレーシングやボリューム散乱シミュレーションのような高価な技術の代わりに、映画「The Matrix Reloaded」の人間の顔のレンダリングのための表面下散乱を近似するために、集められた放射照度テクスチャをブラーするアーティストフレンドリーなカーネルを使用しています。NVIDIA human demo [d’Eon et al., 2007] における d’Eon と彼の同僚による研究は、ダイポール拡散プロファイルのリアルな部分と、ダイポールを近似するためにガウスフィルターを使用するリアルタイムな部分とを、テクスチャー空間に一緒に持ってきました。彼らは、オフラインレンダラーにマッチした、人間の皮膚のための最初のリアルでリアルタイムのサブサーフェイスを作成しました。
拡散プロファイルは、ゼロ平均の等方性ガウス分布の和で近似され、次のようになります。
\begin{eqnarray}
R_d(r) \approx A_g(r) = \sum_{i=1}^k w_i G(v_i, r) \tag{2.48}
\end{eqnarray}
我々は誤差関数を最小化する必要があります [d’Eon et al., 2007] 。
\begin{eqnarray}
V(v_1, \cdots, v_k | R_d) = \int_0^{\infty} r(R_d(r) – \sum_{i=1}^k w_i G(v_i, r))^2 dr \tag{2.49}
\end{eqnarray}
ここで、\(A_g(r)\)はガウス近似、等方性ガウスは\(G(v, r) = \frac{1}{2\pi v}e^{-\frac{r^2}{2v}}\)を持ちます。\(k\)はガウスの数、\(w_i\)は\(\sum_i^k w_i = \int_0^{\infty} 2\pi r R_d(r)dr\)による正規化重みです。式 2.49 を最小化する Levenberg-Marquardt 最適化では，相対 RMS powe 誤差（\(\frac{ \sqrt{V(v, \cdots, v_k | R_d}) }{ \sqrt{\int_0^{\infty} r(R(r))^2 dr } }\) と定義）は \(k=4\) で大理石の緑色波長に対して 1.25%，\(k=8\) で 0.093% の誤差を報告しています [Jensen et al, 2001]。異なる拡散プロファイルの場合、別のガウシアンセットをフィットさせる必要があります。
また、著者は\(k=4\)のときに収束を気にすることなく、より簡単にダイポールにフィットできる閉じた式を発見しました。なぜなら、ダイポールモデルは2つ以上の極関数と分散の和だからです。分散 \(v\) はほぼ一致する。詳細は、原著論文を参照してください。
対称な2次元ガウシアンは1次元の畳み込みに分解されるため、縦と横に分解することができる。GPUでの実装が効率的になります。さらに、複数の放射状プロファイルがある場合、それらを1つずつ新しいガウシアンに結合してから、次のように積分することができます。
\begin{eqnarray}
\sum_{i=1}^{k_1} w_i G(v_i, r) \ast \sum_{j=1}^{k_2} w’_i G(v’_i, r) = \sum_{i=1}^{k_1} \sum_{j=1}^{k_2} w_i w’_j [G(v_i, r) \ast G(v’_j, r)] \tag{2.50}
\end{eqnarray}
2.2.2.2 Pre-integration
表面下散乱は照明が変化するところで最も目立つので，[Penner and Borshukov, 2011]は，オンラインガウシアンサンプリングを事前積分したテクスチャ（照明 のBRDFと影の半影）にオフロードすることによって， 肌のリアルタイム表面下散乱を加速する別の方向に向かっています。シェーディング中は、曲率と影のパラメータに基づいたテクスチャのルックアップを実行するだけで十分です。
スムースサーフェイス. 表面下散乱は、勾配変化領域で最もよく観察されます。法線と影の半影を使用することで、高速に近似することができます。異なる曲率に対してあらかじめ積分された2Dテクスチャ\(T(\frac{1}{r}, N \cdot L)\)は、次のようにフィットさせることができます。
\begin{eqnarray}
D(\theta, r) = \frac{ \int_{-\pi}^{\pi} {\rm max}(0, \cos(\theta + x)) \cdot R(2r \sin(x/2)) dx }{ \int_{-\pi}^{\pi} R(2 \sin(x/2)) dx } \tag{2.51}
\end{eqnarray}
ここで、\(r\)は半径、平面曲率\(\kappa = \frac{1}{r}\)、\(\theta\)は法線\(N\)からのオフセット角度、\(R\)はガウス和としての拡散プロファイルです。このモデルは、マテリアルが滑らかな球面上にあると仮定しています。
細かい詳細. しかし、肌の小さなシワや毛穴のような細かいディテールは、滑らかなメッシュにディテール法線マップで表現されることが多いです。それらの領域に対して正しい結果を近似的に得るために、法線は各チャンネル(R, G, B)の拡散プロファイルで事前にフィルタリングすることができます。4つの法線マップはリアルタイムレンダリングで帯域幅を必要とするため、赤色法線のみが作成され、他の2つはスペキュラー法線とクランプされブラーされた赤色法線の間のブレンドとして作成されています。
影. 影領域については、半影\(T(\frac{1}{w}, s)\)に基づいて事前に積分された2Dテクスチャは、以下のように適合させることができます。
\begin{eqnarray}
P_s(s, w) = \frac{ \int_{\infty}^{\infty} P'(P^{-1}(s) + x) \cdot R(x/w) dx }{ \int_{-\infty}^{\infty} R(x/w) dx } \tag{2.52}
\end{eqnarray}
ここで、\(P(\cdot)\)は影の半影で、光が完全に明るい状態から暗い状態に移行する方法を決定する1次元フォールオフ関数です。これは、半影の幅（完全な光までの距離）をワールド空間における元の影の値\(s\)にマッピングします。\(P^{-1}(\cdot)\) は、半影の幅を問い合わせるための逆関数です。\(P'(\cdot)\)は新しいシフトされた影の半影の関数で、影によるフォールオフの前に表面下散乱の完全な積分を可能にします。このシフトは、表面下散乱を事前積分した2Dテクスチャに保存することができる半影のセグメントを残します。
2.2.2.3 Separable Approximation
低ランク近似. 拡散プロファイルのガウス近似（例えば4つのガウスの和）から個々のガウスを分離し、一つずつぼかしをかけるのではなく［d’Eon et al., 2007］，Jimenezら［2015］はプロファイル全体を特異値分解（SVD）で低ランク近似することを以下のように提案しています。
\begin{eqnarray}
A_s &=& U \sum V^T, \tag{2.53} \\
U &=& (u^{(1)}| u^{(2)} | \cdots | u^{(m)}), \tag{2.54} \\
V &=& (v^{(1)}| v^{(2)} | \cdots | v^{(m)}), \tag{2.55} \\
\sum &=& diag (\sigma_1, \sigma_2, \cdots , \sigma_m), \tag{2.56}
\end{eqnarray}
ここで、\(A_d \in R^{m \times m}\)は極座標における連続プロファイル \(R\) (または便利な実時間近似として \(A_g(r)\)) の直交座標における離散版です。\(u^{(i)}\) と \(v^{(i)}\) は行列の第 \(i\) 列で, \(\sigma\) は降順とします。そして，Eckart-Young の定理 [Eckart and Young, 1936, Stewart, 1993] に基づく \(r\)-rank 近似は，最初の \(r\) 個の最大特異値を残し，その他の値を 0 で置き換えることです．
\begin{eqnarray}
A^{\{r\}}_s &=& U \sum^{\{r\}} V^T \tag{2.57} \\
\sum^{\{r\}} &=& diag (\sigma_1, \cdots, \sigma_r, 0, \cdots, 0) \tag{2.58}
\end{eqnarray}
1次元フィルタを用いてリアルタイムレンダリングのキャッシュコヒーレンス性を高めるために、離散プロファイル\(A^{\{r\}}_d\)はさらに、以下のような可分モデル\(a_i\)に分解されます。
\begin{eqnarray}
A^{\{r\}}_s = \sum_{i=1}^r \sqrt{\sigma_i} u^{(i)} \otimes \sqrt{\sigma_i}v^{(i)} = \sum_{i=1}^r a_i \otimes a_i \tag{2.59}
\end{eqnarray}
ここで、\(\otimes\)は外積です。プロファイル\(R\)は半径方向に対称であるため、\(u^{(i)} = v^{(i)}\)となります。それは1次元フィルター関数\(a_i\)となります。2次元の r-rank フィルタは、水平方向と垂直方向に \(r\) 個の分離可能なフィルタを持つ放射輝度テクスチャに適用することができます。\(i\)番目のパスでは、水平と垂直の両方が同じ重み\(a_i\)を持ちます。
事前積分済みモデル. しかし、放射照度が \(\frac{ \partial b}{ \partial x \partial y} = \frac{ \partial b}{ \partial y \partial x} = 0\) で加法可分 \(b(x,y) = b(x) + b(y)\) であることが分かっている場合、「ランク-1近似」事前積分可分カーネルは
\begin{eqnarray}
A_p(x, y) = \frac{1}{|| R||_1} a_p(x) a_p(y) \tag{2.60}
\end{eqnarray}
であり，式2.38の直交座標で次のように導かれる。
\begin{eqnarray}
L_o(x, y) &=& \int R(x – x’, y – y’) b(x – x’, y – y’) dx’ dy’ \tag{2.61} \\
&=& \int a_p(x – x’) b(x – x’) dx’ + \int a_p(y – y’)b(y – y’) dy’ \tag{2.62} \\
&=& \int \int \frac{1}{||R||_1} a(x – x’) a(y – y’) b(x -x’, y-y’) dx’ dy’ \tag{2.63} \\
&=& \int \int A_p (x – x’, y – y’) b(x – x’, y – y’) dx’ dy’ \tag{2.64} \\
a_p(x) &=& \int R(x, y) dy \tag{2.65} \\
a_p(y) &=& \int R(x, y) dx \tag{2.66} \\
||R||_1 &=& \int a(x) dx = \int a(y) dy \tag{2.67}
\end{eqnarray}
この仮定により、影が垂直または水平のどちらであっても正確な結果を得ることができます。しかし、この定式化は、\(R=R_d\)のため、近似時に固定されている測定されたダイポールモデルのガウス近似であることが多い拡散プロファイルのため、元の論文[Jimenez et al., 2015]で述べたようにアーティストフレンドリーではありません。しかし、リアルタイムレンダリングアプリケーションでは、結果が中心ピクセルと事前積分との間で補間される事前積分可分モデル実装が広く採用されています。Burleyの正規化プロファイルの拡散反射率プロファイルとして代替モデル\(R\)を持っているので，拡散平均自由経路のパラメータ化\(\Theta_3\)を選択すると制御がアーティストフレンドリーになり，この可分モデルと我々の提案アルゴリズムのハイブリッドでより良い品質と性能を実現できます（5.3.2項参照）。