Improved moment shadow maps for translucent occluders, soft shadows and single scattering

拙訳は抜粋や意訳を多く含みます。原著を必ず確認してください。

拙訳

1. はじめに（Introduction）

フィルタ可能なシャドウマップには、半透明な遮蔽物をフィルタ可能なシャドウマップに直接レンダリングする手法[Delalandre et al. 2011Delalandre, C., Gautron, P., Marvie, J.-E. and François, G. 2011. Transmittance function mapping. Symposium on interactive 3D graphics and games 31–38. 10.1145/1944745.1944751. http://gautron.pascal.free.fr/publications/I3D2011/transmittanceFunctionMapping.pdf.; McGuire and Mara 2016McGuire, M. and Mara, M. 2016. A phenomenological scattering model for order-independent transparency. Proceedings of the 20th ACM SIGGRAPH symposium on interactive 3D graphics and games 149–158. 10.1145/2856400.2856418. https://research.nvidia.com/sites/default/files/pubs/2016-02_A-Phenomenological-Scattering/McGuire2016Transparency.pdf.]、PCSS[Fernando 2005Fernando, R. 2005. Percentage-closer soft shadows. ACM SIGGRAPH 2005 sketches 35–es. 10.1145/1187112.1187153. https://developer.download.nvidia.com/shaderlibrary/docs/shadow_PCSS.pdf.]、Summed-Area Table[Crow 1984Crow, F. C. 1984. Summed-area tables for texture mapping. Proceedings of the 11th annual conference on computer graphics and interactive techniques 207–212. 10.1145/800031.808600.]を用いてソフトシャドウを固定コストで行う方法[Lauritzen 2007Lauritzen, A. 2007. Summed-Area Variance Shadow Maps. GPU Gems 3. https://developer.nvidia.com/gpugems/gpugems3/part-ii-light-and-shadows/chapter-8-summed-area-variance-shadow-maps.; Annen et al. 2008Annen, T., Dong, Z., Mertens, T., Bekaert, P., Seidel, H.-P. and Kautz, J. 2008. Real-time, all-frequency shadows in dynamic scenes. ACM Trans. Graph. 27, 3, 1–8. 10.1145/1360612.1360633. https://doi.org/10.1145/1360612.1360633.; Yang et al. 2010Yang, B., Dong, Z., Feng, J., Seidel, H.-P. and Kautz, J. 2010. Variance soft shadow mapping. Computer Graphics Forum 29, 7, 2127–2134. 10.1111/j.1467-8659.2010.01800.x. https://jankautz.com/publications/VSSM_PG2010.pdf.; Shen et al. 2013Shen, L., Feng, J. and Yang, B. 2013. Exponential soft shadow mapping. Computer Graphics Forum 32, 4, 107–116. 10.1111/cgf.12156. http://www.cad.zju.edu.cn/home/jqfeng/papers/Exponential%20Soft%20Shadow%20Mapping.pdf.]、視線に沿った単一散乱の累積をフィルタ可能なシャドウマップに事前計算する手法[Klehm et al. 2014Klehm, O., Seidel, H.-P. and Eisemann, E. 2014. Filter-based real-time single scattering using rectified shadow maps. Journal of Computer Graphics Techniques (JCGT) 3, 3, 7–34. http://jcgt.org/published/0003/03/02/.]、といった多種多様な応用が存在する。

本論ではこれらすべてのアプローチがMoment Shadow Mappingと互換性があることを実証する。その導出の過程で、モーメントベースのブロッカーサーチや6モーメントのMoment Shadow Mappingといった、新しい構成要素を開発する。出来上がるテクニックは魅力的なクオリティ対ランタイムのトレードオフを提供する。また、フラグメントあたりのオーバーヘッドが小さく固定であるため、高出力解像度に合わせて特に良くスケールする。

2. 改善されたMoment Shadow Maps（Improved Moment Shadow Maps）

2.1. 深度分布（Depth Distributions）

フィルタカーネルが $n \in \mathbb{N}$ 個の重み $w_0, ... , w_{n - 1} > 0$ を持ち、フィルタ領域内でサンプルされる対応する深度が $z_0, ... , z_{n-1} \in \mathbb{R}$ であるとすると、これらの量は以下の深度分布を定義する。

Z := \sum_{l = 0}^{n - 1} w_l \cdot \delta_{z_l}

この表記は、深度 $z_l$ を描く確率がフィルタ重み $w_l$ によって求まるようにフィルタ領域を無作為にサンプルする、という確率的な解釈を提案する。そして、各深度 $z$ をあるベクトル量 $\boldsymbol{b}(z) \in \mathbb{R}^{m + 1}$ (ここで、 $m \in \mathbb{N}$ )にマップすることを考える。無作為な深度を入力として用いると、出力の期待値は以下で求まる。

b := \varepsilon_Z(\boldsymbol{b}) := \sum_{l = 0}^{n - 1} w_l \cdot \boldsymbol{b}(z_l) \in \mathbb{R}^{m + 1}

フィルタ可能なシャドウマップはちょうどこの方法で構築される。各テクセルにはある関数 $\boldsymbol{b} : \mathbb {R} \rightarrow \mathbb{R}^{m + 1}$ による $\boldsymbol{b}(z)$ を格納して、上記のフィルタカーネルでシャドウマップをフィルタリングすることで $b = \varepsilon_Z(\boldsymbol{b})$ を生成する。関数 $\boldsymbol{b}$ には $b$ の情報《knowledge》が $Z$ のおおよその再構築を可能にするものが選ばれる。

2.2. 関連研究（Related Work）

Percentage-Closer Filtering[Reeves et al. 1987Reeves, W. T., Salesin, D. H. and Cook, R. L. 1987. Rendering antialiased shadows with depth maps. Proceedings of the 14th annual conference on computer graphics and interactive techniques 283–291. 10.1145/37401.37435. https://doi.org/10.1145/37401.37435.]はシャドウマップから $n \in \mathbb{N}$ 個のサンプルを取ることで｜無理やり《through brute force》深度分布を再構築する。

Z(\boldsymbol{z} < z_f) := \sum_{l = 0}^{n - 1} w_l \cdot \begin{cases} 1 & \text{if } z_l < z_f \\ 0 & \text{otherwise} \end{cases} \text{ where } \boldsymbol{z}(z) := z

Variance Shadow Maps[Donnelly and Lauritzen 2006Donnelly, W. and Lauritzen, A. 2006. Variance shadow maps. Proceedings of the 2006 symposium on interactive 3D graphics and games 161–165. 10.1145/1111411.1111440. https://pierremezieres.github.io/site-co-master/references/vsm_paper.pdf.]は $\boldsymbol{b}(z) := (1, z, z^2)^\mathsf{T}$ を用いる。1つ目の要素は $\varepsilon_Z(\boldsymbol{b_0}) = 1$ なので格納する必要がなく、残りの2つで $Z$ の平均と分散を計算して、Cantelliの不等式で再構築できる。

\mu := b_1, \sigma^2 := b_2 - b_1^2, Z(\boldsymbol{z} < z_f) \ge 1 - \frac{\sigma^2}{\sigma^2 + (z_f - \mu)^2} \text{ if } z_f \ge \mu

同様に、Exponential Shadow Maps[Salvi 2008Salvi, M. 2008. Rendering Filtered Shadows with Exponential Shadow Maps. ShaderX6. http://www.shaderx6.com/index.html.; Annen et al. 2008Annen, T., Mertens, T., Seidel, H.-P., Flerackers, E. and Kautz, J. 2008. Exponential shadow maps. Proceedings of graphics interface 2008 155–161. https://jankautz.com/publications/esm_gi08.pdf.]は $\boldsymbol{b}(z) := (1, \exp(c_{esm} \cdot z))^\mathsf{T}$ (ここで、 $c_{esm} \gg 1$ )を用いて、Markovの不等式で再構築する。

Z(\boldsymbol{z} < z_f) \ge 1 - \frac{b_1}{\exp(c_{esm} \cdot z_f)}

Convolution Shadow Maps[Annen et al. 2007Annen, T., Mertens, T., Bekaert, P., Seidel, H.-P. and Kautz, J. 2007. Convolution shadow maps. Rendering techniques. 10.2312/EGWR/EGSR07/051-060. https://jankautz.com/publications/csmEGSR07.pdf.]は $\boldsymbol{b}$ の関数にフーリエ基底関数を設定し、打ち切られたフーリエ級数を用いる。Exponential Variance Shadow Maps[Lauritzen and McCool 2008Lauritzen, A. and McCool, M. 2008. Layered variance shadow maps. Proceedings of graphics interface 2008 139–146.]はパラメータ $c_{evsm}^+, c_{evsm}^- > 0$ を固定し、2つの指数関数的に｜曲げられた《warped》深度にVariance Shadow Mappingを適用する。すなわち、以下を用いる。

\boldsymbol{b}(z) = (1, \exp(c_{evsm}^+ \cdot z), \exp(c_{evsm}^+ \cdot z)^2, -\exp(c_{evsm}^- \cdot z), \exp(c_{evsm}^- \cdot z)^2)^\mathsf{T}

2.2.1. Moment Shadow Mapping

Moment Shadow Mapping[Peters and Klein 2015Peters, C. and Klein, R. 2015. Moment shadow mapping. Proceedings of the 19th symposium on interactive 3D graphics and games 7–14. 10.1145/2699276.2699277. https://momentsingraphics.de/I3D2015.html.]は以下を用いる。

\boldsymbol{b}(z) = (z^j)_ {j = 0}^m = (1, z, z^2, z^3, z^4)^\mathsf{T}

モーメント数 $m \in \mathbb{N}$ は任意の偶数を取ることができる。

深度 $z_f \in \mathbb{R}$ のフラグメントをシェーディングするとき、シャドウ強度 $Z(\boldsymbol{z} < z_f)$ を推定する必要があるが、シャドウマップからフィルタされたサンプル $b = \varepsilon_Z(\boldsymbol{b})$ のみが与えられている。一般に、これらのモーメントに互換性のある $\mathbb{R}$ 上の深度分布 $S$ は無数に存在する。誤ったセルフシャドウ(サーフェスアクネ)を減らすため、我々は以下で求まる最も小さなシャドウ強度を産む再構築を用いる。

\inf \{ S(\boldsymbol{z} < z_f) | S\text{は}\varepsilon_S(\boldsymbol{b}) = b\text{を持つ}\mathbb{R}\text{上の深度分布} \} \tag{1}

この最小化問題は非自明であるが、すでに解かれている[Kreĭn and Nudel'man 1977Kreĭn, M. G. and Nudel'man, A. A. 1977. The markov moment problem and extremal problems. https://bookstore.ams.org/MMONO/50.]。最適な深度分布 $S$ は深度値 $z_f$ と他の $\frac{m}{2} = 2$ つの深度値のみを用いる。そのような深度分布は唯一であり、アルゴリズム1で効率良く計算する。

アルゴリズム1では、行列 $B(b)$ は正定値《positive definite》でないと正しく動作しない可能性があるが、有意な入力 $b = \varepsilon_Z(\boldsymbol{b})$ に対して半正定値《positive semi-definite》となることが知られている。丸め誤差が入力を不正にする可能性があるが、これらは小さな補間重み $0 < \alpha_b \ll 1$ を用いた固定の定数ベクトルに向かう線形補間によって相殺される。このバイアスは、近似の正確さを低下させたり、ライトリークを増やしたりするので、最小限に留めるべきである。

アルゴリズム1 --- 式(1)を解く
入力: モーメント $b \in \mathbb{R}^{m+1}$ 、フラグメントの深度 $z_f \in \mathbb{R}$
出力: 式(1)での表現を最小化する深度分布 $S$
1. $B(b) := (b_{j+k})_ {j,k=0}^{\frac{m}{2}} = \left( \begin{array}{c} b_0 & b_1 & \cdots & b_{\frac{m}{2}} \\ b_1 & b_2 & \cdots & b_{\frac{m}{2}+1} \\ \vdots & \vdots & \ddots & \vdots \\ b_{\frac{m}{2}} & b_{\frac{m}{2}+1} & \cdots & b_m \end{array} \right) \in \mathbb{R}^{(\frac{m}{2}+1)\times(\frac{m}{2}+1)}$ をセットする。
2. $q \in \mathbb{R}^{\frac{m}{2}+1}$ に対する $B(b) \cdot q = (1, z_f^1, \dots, z_f^{\frac{m}{2}})^\mathsf{T}$ を解く。
3. $z$ に対する多項式 $\sum_{j=0}^{\frac{m}{2}} q_j \cdot z^j = 0$ を解き、 $z_1, \dots , z_{\frac{m}{2}} \in \mathbb{R}$ による個別の解を示す。
4. $z_0 := z_f$ と $A := (z_l^j)_ {j,l=0}^{\frac{m}{2}} = \left( \begin{array}{c} z_0^0 & \cdots & z_{\frac{m}{2}}^0 \\ \vdots & \ddots & \vdots \\ z_0^{\frac{m}{2}} & \cdots & z_{\frac{m}{2}}^{\frac{m}{2}} \end{array} \right) \in \mathbb{R}^{(\frac{m}{2}+1)\times(\frac{m}{2}+1)}$ をセットする。
5. $w \in \mathbb{R}^{\frac{m}{2}+1}$ に対する $A \cdot w = (b_0, b_1, \dots, b_{\frac{m}{2}})^\mathsf{T}$ を解く。
6. $\sum_{l=0}^{\frac{m}{2}} w_l \cdot \delta_{z_l}$ を返す。

4つのモーメントにそれぞれ16ビットを単に充てるとすると丸め誤差が強く出てしまうので、出力データ型の表現できる範囲を侵さないで行列式を最大化する数値的な最適化となるようなアフィン変換を施してから16ビット浮動小数点数に格納する。

2.3. 符号付き深度（Signed Depth）

元々[Peters and Klein 2015Peters, C. and Klein, R. 2015. Moment shadow mapping. Proceedings of the 19th symposium on interactive 3D graphics and games 7–14. 10.1145/2699276.2699277. https://momentsingraphics.de/I3D2015.html.]では深度を区間 $[0, 1]$ で定義するよう提案したが、この選択は最適な選択である $[-1, 1]$ より悪いことが判明した。このときニア面が $-1$ 、ファー面が $1$ に位置する。

この定義を支持する理由はモーメントでの線形変換の効果に基づいている。 $x \ne 0$ である $x, y \in \mathbb{R}$ が深度値の線形変換 $x \cdot z + y$ を記述すると、以下のモーメントに対する線形変換を導く。

\varepsilon_Z((x \cdot \boldsymbol{z} + y)^j) = \varepsilon_Z \left( \sum_{k=0}^j \left( \begin{array}{c} j \\ k \end{array} \right) \cdot (x \cdot \boldsymbol{z})^k \cdot y^{j - k} \right) = \sum_{k=0}^j \left( \begin{array}{c} j \\ k \end{array} \right) \cdot x^k \cdot y^{j-k} \cdot b_k \tag{2}

この線形変換は対角成分 $x^0, \dots, x^m$ を持つ下三角行列に対応する。すなわち、この行列式は $\prod_{j=1}^m x^j$ である。

区間 $[-1, 1]$ で深度を定義すると仮定するとして、上記の変換を適用すると、 $j$ 番目のモーメントは最大で $(|x|+|y|)^j$ の大きさを取り得る。我々は浮動小数点数を用いているので、大きさの最大値を1にするためにこの大きさで割る場合には相対的な精度は減少しない。これを先程の変換と組み合わせると、その行列式は以下のようになる。

\prod_{j=1}^m \left( \frac{x}{|x|+|y|} \right)^j

明らかに、この行列式の大きさは $y = 0$ のときに最大であり、 $x$ の変化に対して不変である。

式(2)の変換は特殊な量子化変換のように働く。その行列式を最大化することはMoment Shadow Mapに格納されるデータのエントロピーに良い影響を持つ。故に、深度区間 $[-1, 1]$ の選択が確かに最適である。 $m = 4$ の場合、 $[0, 1]$ から $[-1, 1]$ に変えると行列式は1024になる。これはおよそ $\log_2 1024 = 10$ ビットのエントロピーの増加に対応する。

2.4. 疎な量子化変換（Sparse Quantization Transform）

16ビットの固定小数点チャネルによるMoment Shadow Mapを用いる場合、格納することで発生する丸め誤差はアルゴリズム1に起因する丸め誤差をはるかに上回る。量子化変換を用いることで、この丸め誤差を最高効率で最小化する。故に、符号付き深度を用いても精度は改善しないが、最適化を行うことができるようになる。

オリジナルの量子化変換の適用コストは著しく、これはMoment Shadow Mapのテクセルごとやシェードされたフラグメントごとに密な4x4行列によるベクトルの乗算を必要とするが、 $[-1, 1]$ での深度値に起因するモーメントを直接格納するのと比べて、この変換はエントロピーを4.28ビット増加させる。

我々は変換の最適化における探索空間をエントロピー的に大きく減らすことなく実質的に制限することができることを発見した。特に、以下の変換はエントロピーを4.21ビット増加させる。

\Theta_4^* \left( \begin{array}{c} b_1 \\ b_2 \\ b_3 \\ b_4 \end{array} \right) := \left( \begin{array}{c} \frac{3}{2} & 0 & -2 & 0 \\ 0 & 4 & 0 & -4 \\ \frac{\sqrt{3}}{2} & 0 & -\frac{2\sqrt{3}}{9} & 0 \\ 0 & \frac{1}{2} & 0 & \frac{1}{2} \end{array} \right) \cdot \left( \begin{array}{c} b_1 \\ b_2 \\ b_3 \\ b_4 \end{array} \right) + \left( \begin{array}{c} \frac{1}{2} \\ 0 \\ \frac{1}{2} \\ 0 \end{array} \right) \tag{3}

故に、結果の品質は実践的には変化しないが、成分の半数が0であるので、この変換の適用は大まかに2倍速くなる。

2.5. ワーストケースに対するバイアス（Biasing for the Worst Case）

オリジナルのスキームは以下によってモーメントのベクトル $b \in \mathbb{R}^5$ を置き換える。

b' := (1 - \alpha_b) \cdot b + \alpha_b \cdot b^* \text{ where } b^* := (1, 0, 1, 0, 1)^\mathsf{T} \in \mathbb{R}^5

このスキームは平均的なケースで最良の結果をもたらすために導入された[Peters and Klein 2015Peters, C. and Klein, R. 2015. Moment shadow mapping. Proceedings of the 19th symposium on interactive 3D graphics and games 7–14. 10.1145/2699276.2699277. https://momentsingraphics.de/I3D2015.html.]。

しかし、これが失敗するようなワーストケースが存在する。関連するシャドウキャスタをクリップせずに深度範囲をより活用するために深度値を区間 $[-1, 1]$ にクランプする場合、以下のような深度分布の形式が生じる。

Z = (1 - w_1) \cdot \delta_{-1} + w_1 \cdot \delta_1

$b^* = \varepsilon_{Z^* }(\boldsymbol{b}), Z^* = \frac{1}{2}\delta_{-1} + \frac{1}{2}\delta_{1}$ とすると、バイアスされたベクトル $b'$ は以下の深度分布に対応する。

(1 - \alpha_b) Z + \alpha_b Z^* = \left( (1 - \alpha_b)(1 - w_1) + \frac{\alpha_b}{2} \right)\delta_{-1} + \left( (1 - \alpha_b)w_1 + \frac{\alpha_b}{2} \right) \delta_1

明らかに、この深度分布は $-1$ と $1$ の2つの深度値のみを用いる。そのような深度分布は有効な領域の境界でのモーメントのベクトルに対応する。すると、多少の丸め誤差でも無効な領域に入ってしまい、バイアス処理はその目的を達成できなくなる。つまり、深度値をクランプすると、ある場面で不安定性が現れるようになる。

すべてのケースでロバストに振る舞うバイアスを導入するため、ワーストケースにおける最大効率、すなわち、有効な領域 $\text{conv } \boldsymbol{b}(\mathbb{R})$ の境界への最大の距離を持つモーメントのベクトル $b^* \in \text{conv } \boldsymbol{b}([-1, 1])$ を求める。量子化変換を用いない場合、最適化結果は以下になる。

b^* = (1, 0, 0.375, 0, 0.375)^\mathsf{T}

モーメントを単精度の浮動小数点数で格納すると、モーメントバイアスが $\alpha_b = 3\cdot10^{-7}$ あれば十分であることを発見した。

式(3)の量子化変換を用いるとき、変換されたモーメントに強い丸め誤差が載るため、 $\Theta_4^*$ を適用してから計算を行うべきである。この場合、最適化結果は以下になる。

b^* = (1, 0, 0.628, 0, 0.628)^\mathsf{T}

このとき、テクセルあたり64ビットを用いると、モーメントバイアスが $\alpha_b = 6 \cdot 10^{-5}$ あれば確実にアーティファクトを排除することを観測した。

2.6. 結果と議論（Results and Discussion）

3. 半透明遮蔽物のシャドウ（Shadow for Translucent Occluders）

3.1. 関連研究（Related Work）

半透明な遮蔽物が存在するとき、ライトから表面への透過率は複雑な方法で深度に依存することがある。これは区分線形関数で表されるが、そのフィルタリングは自明ではない[Salvi et al. 2010Salvi, M., Vidimče, K., Lauritzen, A. and Lefohn, A. 2010. Adaptive volumetric shadow maps. Computer Graphics Forum 29, 4, 1289–1296. 10.1111/j.1467-8659.2010.01724.x. https://www.intel.com/content/dam/develop/external/us/en/documents/salvi-avsm-egsr2010-564567.pdf.]。他のアプローチではPercentage-Closer Filteringを用いたり、半透明率に比例してフラグメントを無作為に破棄したりする[Enderton et al. 2010Enderton, E., Sintorn, E., Shirley, P. and Luebke, D. 2010. Stochastic Transparency. Proceedings of the 2010 ACM SIGGRAPH symposium on interactive 3D graphics and games 157–164. 10.1145/1730804.1730830. https://luebke.us/publications/StochasticTransparency_I3D2010.pdf.; McGuire and Enderton 2011McGuire, M. and Enderton, E. 2011. Colored stochastic shadow maps. Symposium on interactive 3D graphics and games 89–96. 10.1145/1944745.1944760. https://casual-effects.com/research/McGuire2011CSSM/index.html.]。

Fourier Opacity Mapping[Jansen and Bavoil 2010Jansen, J. and Bavoil, L. 2010. Fourier opacity mapping. Proceedings of the 2010 ACM SIGGRAPH symposium on interactive 3D graphics and games 165–172. 10.1145/1730804.1730831. https://volumetricshadows.wordpress.com/wp-content/uploads/2011/06/fourier-opacity-mapping.pdf.]は、Convolution Shadow Maps[Annen et al. 2007Annen, T., Mertens, T., Bekaert, P., Seidel, H.-P. and Kautz, J. 2007. Convolution shadow maps. Rendering techniques. 10.2312/EGWR/EGSR07/051-060. https://jankautz.com/publications/csmEGSR07.pdf.]のアイデアを導入して、ソートせずに｜足し合わせ《be accumulated additively》られる吸収関数(透過率の対数)をフーリエ級数で表現する。Translucent Shadow Maps[Delalandre et al. 2011Delalandre, C., Gautron, P., Marvie, J.-E. and François, G. 2011. Transmittance function mapping. Symposium on interactive 3D graphics and games 31–38. 10.1145/1944745.1944751. http://gautron.pascal.free.fr/publications/I3D2011/transmittanceFunctionMapping.pdf.]は似たようなアプローチを取るが、ソートを必要とする透過率関数を表現する。Phenomeonlogical Scattering[McGuire and Mara 2016McGuire, M. and Mara, M. 2016. A phenomenological scattering model for order-independent transparency. Proceedings of the 20th ACM SIGGRAPH symposium on interactive 3D graphics and games 149–158. 10.1145/2856400.2856418. https://research.nvidia.com/sites/default/files/pubs/2016-02_A-Phenomenological-Scattering/McGuire2016Transparency.pdf.]はVariance Shadow Mapを通じて透過率を表現するが、フラグメントを確率的に破棄することでソートを回避する。このテクニックはヒューリスティックなコースティクスも追加する。

3.2. 半透明遮蔽物に対するMoment Shadow Maps（Moment Shadow Maps for Translucent Occluders）

我々のアプローチはMoment Shadow Mapが透過率を表現すると言う点でTranslucent Shadow Mapsに似ている。吸収関数の表現には全吸収《total absorption》 $b_0$ のための追加のチャネルを必要とするだろう。さらに、不透明遮蔽物が無限の吸収に対応するとして、不透明と半透明の遮蔽物に対する単一のMoment Shadow Mapを用いたいと考えている。

この選択の欠点はMoment Shadow MapにレンダリングするときにOrder-Independent Transparencyのための手法を必要とすることである。我々はこれが我々の貢献と無関係《orthogonal》であると考えており、既存の手法(例えば、Stochastic Transparency[Enderton et al. 2010Enderton, E., Sintorn, E., Shirley, P. and Luebke, D. 2010. Stochastic Transparency. Proceedings of the 2010 ACM SIGGRAPH symposium on interactive 3D graphics and games 157–164. 10.1145/1730804.1730830. https://luebke.us/publications/StochasticTransparency_I3D2010.pdf.])は上手く動作するはずである。実験ではソートされたジオメトリに依存している。

不透明度 $\alpha_0, \dots , \alpha_{n_s - 1} \in [0, 1]$ を持つ深度 $z_0 < z_1 < \dots < z_{n_s - 1}$ における光線に沿った $n \in \mathbb{N}$ 個の表面があるとする。深度 $z_f \in \mathbb{R}$ へと透過した光量は関連する透過率係数の積 $\prod_{k=0,z_k<z_f}^{n_s-1} (1-\alpha_k)$ である。この透過率は、深度 $z_j$ において、残りの光の割合 $\alpha_j$ がブロックされるので、以下の深度分布によって正確にモデル化される。

Z := \sum_{j=0}^{n_s-1} \left( \prod_{k=0}^{j-1} 1-\alpha_k \right) \alpha_j \delta_{z_j}

ブレンディングのover operator¹を用いてMoment Shadow Mapに後から前へこれらの表面をレンダリングすると仮定する。Moment Shadow Mapは適宜クリアするので、 $z_{n_s - 1} = \alpha_{n_s - 1} = 1$ と仮定するのは安全である。最終的に、モーメントのベクトルは以下のようになる。

b := \sum_{j=0}^{n_s-1} \left( \prod_{k=0}^{j-1} 1 - \alpha_k \right) alpha_j \boldsymbol{b}(z_j) = \varepsilon_Z(\boldsymbol{b})

近似の誤差はアルゴリズム1を通してモーメントからシャドウ強度を再構築する間にのみ引き起こされる。注意点として、 $b_0$ はブレンドし合ったすべての表面のアルファの合計に対応しており、 $\alpha_{n_s - 1} = 1$ であることから $b_0 = 1$ であると分かっているので、依然として $b_0$ は格納する必要はない。

over operatorが必要とされるので、半透明遮蔽物はMoment Shadow Mapに直接レンダリングしなければならない。さらに、現代のグラフィクスAPIの制約から、ブレンディングで用いる不透明度はRGBAテクスチャに書き込まれる値から独立することはできないので、代わりに各チャネルが16ビットのRGテクスチャを2つ用いる。レンダリングはMRTのハードウェアサポートを用いるため、パフォーマンス低下は軽度である。もちろん、疎な量子化変換を用いることも有効である。

3.3. 結果と議論（Results and Discussion）

4. Moment Soft Shadow Mapping

4.1. 関連研究（Related Work）

我々のアプローチはPercentage-Closer Soft Shadows[Fernando 2005Fernando, R. 2005. Percentage-closer soft shadows. ACM SIGGRAPH 2005 sketches 35–es. 10.1145/1187112.1187153. https://developer.download.nvidia.com/shaderlibrary/docs/shadow_PCSS.pdf.]のフレームワークを用いる。PCSSは多少の目立つアーティファクトを含みつつもっともらしい結果を生成するが、過度なサンプリングによりそのコストは高い。

Summed-Area Variance Shadow Maps[Lauritzen 2007Lauritzen, A. 2007. Summed-Area Variance Shadow Maps. GPU Gems 3. https://developer.nvidia.com/gpugems/gpugems3/part-ii-light-and-shadows/chapter-8-summed-area-variance-shadow-maps.]はSummed-Area Tableを用いてそのサンプリングを回避しようとする。Summed-Area Tableを用いると、任意の矩形上の積分が4つのサンプルのみで行うことができるので、矩形エリアライトに対応するボックスカーネルを用いたフィルタリングをそのフィルタサイズに依らずに定数時間で行うことができるようになる。ただし、ブロッカーサーチには依然としてサンプリングを必要とする。

Variance Soft Shadow Mapping[Yang et al. 2010Yang, B., Dong, Z., Feng, J., Seidel, H.-P. and Kautz, J. 2010. Variance soft shadow mapping. Computer Graphics Forum 29, 7, 2127–2134. 10.1111/j.1467-8659.2010.01800.x. https://jankautz.com/publications/VSSM_PG2010.pdf.]はSummed-Area Variance Shadow Mapへの単一クエリに基づくヒューリスティックを用いてブロッカーサーチを高速化する。難しいシチュエーションはより高価な方法で対処される。Convolution Soft Shadow Mapping[Annen et al. 2008Annen, T., Dong, Z., Mertens, T., Bekaert, P., Seidel, H.-P. and Kautz, J. 2008. Real-time, all-frequency shadows in dynamic scenes. ACM Trans. Graph. 27, 3, 1–8. 10.1145/1360612.1360633. https://doi.org/10.1145/1360612.1360633.]はSummed-Area Tableかミップマップのいずれかを用いて、Convolution Shadow Mapsに基づくフィルタリングを行う。ブロッカーサーチにはフィルタ可能な第二のテクスチャセットを用いる。Exponential Soft Shadow Mapping[Shen et al. 2013Shen, L., Feng, J. and Yang, B. 2013. Exponential soft shadow mapping. Computer Graphics Forum 32, 4, 107–116. 10.1111/cgf.12156. http://www.cad.zju.edu.cn/home/jqfeng/papers/Exponential%20Soft%20Shadow%20Mapping.pdf.]はExponential Shadow Mapの小さな領域上でSummed-Area Tableを用いることで、許容できない精度ロスを回避する。これもブロッカーサーチに追加のテクスチャを必要とする。この著者らはカーネルサブディビジョンを用いてガウシアンフィルタカーネルのより良い近似を行う。

4.2. 4モーメントでのSummed-Area Table（Summed-Area Tables with Four Moments）

我々のテクニックはVariance Soft Shadow Mappingと同じ基本工程に従うが、Summed-Area Tableへのクエリを2つ以上必要としない。

Summed-Area Tableは Klehm et al. [2014Klehm, O., Seidel, H.-P. and Eisemann, E. 2014. Filter-based real-time single scattering using rectified shadow maps. Journal of Computer Graphics Techniques (JCGT) 3, 3, 7–34. http://jcgt.org/published/0003/03/02/.] が推奨するように行/列あたり1スレッドを用いる2パスのコンピュートシェーダによって生成される。

Summed-Area Tableはその精度を十分に取るため、最大カーネルサイズについての先行知識[Lauritzen 2007Lauritzen, A. 2007. Summed-Area Variance Shadow Maps. GPU Gems 3. https://developer.nvidia.com/gpugems/gpugems3/part-ii-light-and-shadows/chapter-8-summed-area-variance-shadow-maps.]を活用して、32ビット整数とモジュラー演算を用いる。

用いられる最大のカーネルが $n_t \in \mathbb{N}$ テクセル(例えば、28x28のカーネルで $n_t = 784$ )をカバーすると仮定する。Moment Shadow Mapに格納された変換済みモーメント $\Theta_4^* (b_1, b_2, b_3, b_4)$ は初めは $[0, 1]^4$ にある。これらに $\frac{2^{32}-1}{n_t}$ を乗算して整数に丸め込む場合、関連する最大のカーネルでのモーメントの総和は $\{0, \dots, 2^{32}-1\}$ にあることが知られている。故に、この数値は32ビットの符号なし整数で表現できる。それと同時に、依然として以下の精度を持っている。

\log_2 \frac{2^{32} - 1}{n_t}

これは上記の例では22.4ビットになる。これは単精度floatより若干劣る程度であり、モーメントバイアスが $\alpha = 6 \cdot 10^{-7}$ あれば十分であることを発見した。カーネルが大きくなれば、この値も大きくする必要がある。

我々の実装では、128ビットのMoment Shadow Mapと整数のSummed-Area Tableを生成する。このステップ中に、度々オーバーフローが発生するだろうが、これらは $2^{32}$ の倍数を減算するだけなので、安全に無視できる。 $n_t$ かそれより少ないテクセルを含むカーネルでSummed-Area Tableをクエリするとき、その結果は $\{0, \dots, 2^{32}-1\}$ になければならないことを知っており、故に、整数で計算しても必然的に正しい結果を導く。

Summed-Area Tableを持つことで、ミップマッピングが必要なくなる。従って、64ビットのMoment Shadow Mappingに比べてメモリ必要量は50%の増加するだけになる。

\frac{128 \text{ bits}}{64 \text{ bits} \cdot \frac{4}{3}} = \frac{6}{4} = 150\%

4.3. ブロッカーサーチ（Blocker Search）

ブロッカーサーチでは、サーチ領域の4つのモーメントをクエリするためにSummed-Area Tableで単一のルックアップを行う。そして、アルゴリズム1を用いて、バイアスされたモーメントとバイアスされていないフラグメント深度 $z_f$ を、確率 $w_0, w_1, w_2 > 0$ を持つ3つの深度値 $z0 = z_f, z_1, z_2 \in \mathbb{R}$ で構成される深度分布 $Z := \sum_{l=0}^2 w_l \cdot \delta_{z_l}$ に合わせる。

我々はこの再構築がground truthと一致することを仮定する。サーチ領域が1つか2つの表面を含む場合、再構築はほぼ完璧であることが知られている[Peters and Klein 2015, Proposition 4Peters, C. and Klein, R. 2015. Moment shadow mapping. Proceedings of the 19th symposium on interactive 3D graphics and games 7–14. 10.1145/2699276.2699277. https://momentsingraphics.de/I3D2015.html.]。サーチ領域が表面を3つ含み、そのひとつが深度 $z_f$ にある場合、モーメントにより一意に定まり、正しく再構築される[Peters and Klein 2015, Proposition 8Peters, C. and Klein, R. 2015. Moment shadow mapping. Proceedings of the 19th symposium on interactive 3D graphics and games 7–14. 10.1145/2699276.2699277. https://momentsingraphics.de/I3D2015.html.]。このため、バイアスされていないフラグメント深度を用いることは有益である。これより複雑なケースは稀である。

この分布は仮定により正しいので、Percentage-Closer Soft Shadowsと同様の平均ブロッカー深度を導くことができる。

\frac{\sum_{l=1, z_l < z_f}^2 w_l \cdot z_l}{\sum_{l=1, z_l < z_f}^2 w_l}

しかし、この数式はロバストではない。この除数はちょうどサーチ領域で計算されるシャドウ強度 $Z(\boldsymbol{z} < z_f)$ であり、いくらでも小さくなる可能性がある。この場合、式が意味をなさなくなる。小さいシャドウ強度は遮蔽されないフラグメントを暗に示し、そのようなフラグメントでは、ブロッカーサーチは小さなフィルタカーネルが使われることを示す $z_0 = z_f$ を返すだろう。

この堅牢性を盛り込むために、平均ブロッカー深度を以下とする。

\frac{\varepsilon_{z_0} \cdot z_0 + \sum_{l=1, z_l < z_f}^2 w_l \cdot z_l}{\varepsilon_{z_0} + \sum_{l=1, z_l < z_f}^2 w_l} \tag{4}

ここで、 $\varepsilon_{z_0}$ は0以上の小さな定数である。このパラメータは、一般的に｜完全に照らされた《fully lit》フラグメントに影響を与えるためにその最終結果が変化しないので、品質とって重要ではないことが判明した。我々はすべての実験で $\varepsilon_{z_0} = 10^{-3}$ を用いる。

4.4. フィルタリング（Filtering）

ブロッカー深度が分かれば、半影推定[Fernando 2005Fernando, R. 2005. Percentage-closer soft shadows. ACM SIGGRAPH 2005 sketches 35–es. 10.1145/1187112.1187153. https://developer.download.nvidia.com/shaderlibrary/docs/shadow_PCSS.pdf.]は十分なフィルタサイズを提供する。これとシャドウマップのフラグメントのテクスチャ座標を組み合わせると、フィルタ領域の左上と右下のテクスチャ座標を計算できる。これは一般にテクセルの中心にないため、整数のSummed-Area Tableに対する補間処理が必要になる。

大方の予想とは裏腹に、隣接テクセルの基礎となる値は不明の $2^{32}$ の倍数によって異なるので、フィルタ領域の4つ角のサンプルを直接バイリニア補間するのは正しくない。このような問題は $n_t$ より少ないテクセルを含む領域上の積分で例外的に処理することで回避できる。すなわち、安全に変換できる9つの領域に分割して、その結果をフィルタ領域との交差率で重み付けして総和する。この処理は確実に動作するが、フラグメントあたり $4 \cdot 4 \cdot 4 \cdot 32 = 2048$ ビットを読み込む必要があるので、そのコストはかなり大きい。代替案として、無作為化したディザリングを試したが、ハードシャドウの境界でノイズが強くなりすぎることが判明した。

フィルタサイズは大きくなる可能性があるため、十分な深度バイアスを用いることが重要である。我々はフィルタサイズに比例して増加させることを推奨する。加えて、適応的深度バイアスが用いられるかもしれない[Dou et al. 2014Dou, H., Yan, Y., Kerzner, E., Dai, Z. and Wyman, C. 2014. Adaptive depth bias for shadow maps. Journal of Computer Graphics Techniques (JCGT) 3, 4, 146–162. http://jcgt.org/published/0003/04/08/.]。

4.5. 最適化（Optimization）

テクニックを最適化する最も効率的な方法はテクスチャ読み出しの数を減らすことである。ブロッカーサーチの間の補間コストを回避するため、我々はテクセルグリッドに一致するようにサーチ領域を拡大する。これはソフトシャドウに小さな不連続性を引き起こすが、相当なスピードアップにより容易に正当化される。

他には、フラグメントが本影にあると明らかになったときにフィルタリングをスキップする、という方法がある。シャドウ強度 $\sum_{l=1,z_l<z_f}^2 w_l$ は式(4)の一部としてすでに計算済みであるので、しきい値 $1 - \varepsilon_u$ を上回れば、フラグメントは本影にあるとしてすぐに最大シャドウ強度を返す。我々の実験では、99%以下のシャドウ強度の部分に紐づく $\varepsilon_u = 0.01$ が大きく接続された領域でのテクスチャ読み出しの必要を減少させつつロバストな結果を生み出す。

我々は完全に照らされたフラグメントのフィルタリングステップをスキップしようともしたが、Moment Shadow Mappingでもたらされる下界《lower bound》は多すぎる偽陽性を引き起こす。代わりに上界《upper bound》を用いることもできるが、ほとんどのフラグメントが完全に照らされたとは分類されない。したがって、このアプローチは推奨しないし、実験では用いていない。

4.6. 結果と議論（Results and Discussion）

4.6.1. 定質的評価（Qualitative Evaluation）

4.6.2. 実行時間（Run Time）

4.6.3. 結論（Conclusions）

5. モーメント付き事前フィルタ済み単一散乱（Prefiltered Single Scattering with Moments）

5.1. 関連研究（Related Work）

リアルタイム散乱のアプローチの多くは、そのコストを合理的な範囲に収めるため、大域照明効果を無視し、その代わりに単一散乱にフォーカスする。ピクセルシェーダでのレイマーチング[Toth and Umenhoffer 2009Toth, B. and Umenhoffer, T. 2009. Real-time volumetric lighting in participating media. Eurographics 2009 - short papers. 10.2312/egs.20091048. https://diglib.eg.org/server/api/core/bitstreams/f55f92ba-2775-4dd3-be2b-76a574e50559/content.]はこれを計算する即時的な方法である。レイサンプルごとに、一般的なシャドウマップが可視性を定めるためにサンプルされる。ボクセル化されたシャドウボリューム[Wyman 2011Wyman, C. 2011. Voxelized shadow volumes. Proceedings of the ACM SIGGRAPH symposium on high performance graphics 33–40. 10.1145/2018323.2018329. https://cwyman.org/papers/hpg11_voxelShadows.pdf.]は128クエリを1度に行えるように異なる座標系で二値ボリュームとしてシャドウマップを格納する。1Dのmin-maxミップマップは単一ステップで完全に照らされたか完全に影になったレイセグメントを横断するのに使われるかもしれない[Chen et al. 2011Chen, J., Baran, I., Durand, F. and Jarosz, W. 2011. Real-time volumetric shadows using 1D min-max mipmaps. Symposium on interactive 3D graphics and games 39–46 PAGE@7. 10.1145/1944745.1944752. https://groups.csail.mit.edu/graphics/mmvs/mmvs.pdf.]。

5.1.1. 事前フィルタされた単一散乱（Prefiltered Single Scattering）

事前フィルタされた単一散乱[Klehm et al. 2014Klehm, O., Seidel, H.-P. and Eisemann, E. 2014. Filter-based real-time single scattering using rectified shadow maps. Journal of Computer Graphics Techniques (JCGT) 3, 3, 7–34. http://jcgt.org/published/0003/03/02/.]は単一散乱へのフィルタ可能なシャドウマップの概念を導入する。すべての光線に対する単一散乱の結果はConvolution Shadow Mapに事前計算される。この手法は高速でシーンの複雑さに依存しないが、使われるフーリエ級数がリンギングを引き起こしたり、メモリ必要量が高かったり(例えば、テクセル当たり256ビット)する。

単一散乱の制約に加え、事前フィルタされた単一散乱はいくつかの追加の単純化仮定を盛り込む。関与媒質は均質でなければならない。すなわち、透過率を定義する消散《extinction》係数 $\sigma_t$ 、BRDFのボリューメトリック的類似物である位相関数 $f(\omega_l, \omega_p)$ 、散乱アルベド $\rho$ 、といった物理特性はどこでも同じでなければならない。また、 $\omega_l$ と最大放射照度 $E_l$ を持つ単一のディレクショナルライトからの均質なライティングを仮定する。複数のディレクショナルライトは｜重ね合わせ《superposition》によって扱うことができる。

カメラから距離 $s$ にあり、放射輝度 $L_s$ を出射する表面要素を考える。そのカメラは位置 $p$ にあり、方向 $\omega_p$ を向くとする。3D空間の位置が照らされていれば $1$ 、遮られていれば $0$ に対応する、光の可視性関数を $V(q)$ とすると、カメラが受け取る放射輝度は以下になる。

\exp(-\sigma_t \cdot s) \cdot L_s + f(\omega_l, \omega_p) \cdot E_l \cdot \int_0^s \exp(-\sigma_t \cdot t) \cdot V(p - t \cdot \omega_p) dt

第1の被加数《summand》²は吸収とout-scatteringの後に残る表面からの放射輝度である。第2の被加数は単一散乱をモデル化する。視線に沿う照らされたセグメントごとに、 $f(\omega_l, \omega_p) \cdot E_l$ の差分《differential》放射輝度が足されるが、その $\exp(-\sigma_t \cdot t)$ はカメラに透過される。

単一散乱の計算コストは可視性項を含む積分にある。光線に沿った一次元関数としてみるとき、可視性関数は、最初の遮蔽物までを $1$ 、その以降を $0$ とする、単純なHeavisideの階段関数である。節2.2で述べたフィルタ可能なシャドウマップはフィルタの適用を可能にする方法においてそのような関数を格納する方法をもたらす。我々は単一散乱の積分をシャドウマップの行《rows》の上での積分に変えるのにこれらを活用する。

この終わりに、シャドウマップのパラメータ化は、光線がシャドウマップでの行に対応すること、シャドウマップ内の光線が一定であること、の2つの要件を満足しなければならない。ほとんどの場合、そのようなパラメータ化は単純な透視投影変換として構築できる[Klehm et al. 2014Klehm, O., Seidel, H.-P. and Eisemann, E. 2014. Filter-based real-time single scattering using rectified shadow maps. Journal of Computer Graphics Techniques (JCGT) 3, 3, 7–34. http://jcgt.org/published/0003/03/02/.]。我々はこの線形修正《rectification》を実装したが、節5.2で与えられる理由により、我々は他の提案されたソリューションとして、リサンプリングによって適用される非線形修正変換[Baran et al. 2010Baran, I., Chen, J., Ragan-Kelley, J., Durand, F. and Lehtinen, J. 2010. A hierarchical volumetric shadow algorithm for single scattering. ACM Trans. Graph. 29, 6. 10.1145/1882261.1866200. https://groups.csail.mit.edu/graphics/volumeshadows/vs.pdf.]を選択した。

ワールド空間から修正された空間に座標を変換するため、まずライトビュー空間に変換し、座標系の原点をカメラ位置に移動する。この空間では、ライト方向はZ軸に対応し、他の軸は直交しながらも任意に選択される。シャドウマップの水平テクスチャ座標はXY平面への射影後の原点への距離 $r$ に対応する。これは光線とカメラの間の距離に一致する。他の2つの座標では、球面座標に変換する。垂直テクスチャ座標は方位角《azimuth》 $\varphi \in (0, 2\pi]$ に対応する。シャドウマップに格納された深度は反転した傾斜角《inclination》 $\pi - \theta \in [0, \pi]$ に対応する。この方法では、新しい座標系での深度値はオリジナルの深度値に伴い単調に増加する。

$\varphi$ と $\theta$ はカメラへの距離に独立であるので、視線はシャドウマップの行に対応し、必要な通りに一定の深度を持つ。同時に、このパラメータ化は、各光線が一定の $r$ と $\varphi$ を持つ、すなわち、単一のテクセルに対応するので、シャドウマップで有効である。epipolar幾何の観点では、 $\varphi$ はライト方向を含み、カメラを通過するepipolarスライスを指す。異なるepipolarスライスに対する単一散乱の計算は独立である。シャドウマップにタイトに合わせるため、 $r$ 、 $\varphi$ 、 $\theta$ の境界は視錐台全体をカバーするように計算される(付録A参照)。 $\theta$ 次元に沿って、ライトリークを回避するためにガードバンドを追加する。 $\theta$ の値間隔の長さは10%伸長される。

我々は整数のテクセルインデックス $u, v$ で指されるこの座標系でフィルタ可能なシャドウマップ $b(u, v)$ を生成する。各テクセルは光線に沿った可視性関数の表現(例えば、フーリエ係数[Klehm et al. 2014Klehm, O., Seidel, H.-P. and Eisemann, E. 2014. Filter-based real-time single scattering using rectified shadow maps. Journal of Computer Graphics Techniques (JCGT) 3, 3, 7–34. http://jcgt.org/published/0003/03/02/.]やモーメント)を格納し、行のフィルタリングは対応するepipolarスライスでのすべての視線に沿ったフィルタリングに対応する。関連する積分を事前計算するため、スライスあたりの後続の光線サンプル間のワールド空間の距離 $\Delta(v)$ が分かっている必要がある。この量は傾斜角 $\theta$ に依存するので、ヒューリスティックが必要になる。洗練されたヒューリスティクスが提案されていた[Klehm et al. 2014Klehm, O., Seidel, H.-P. and Eisemann, E. 2014. Filter-based real-time single scattering using rectified shadow maps. Journal of Computer Graphics Techniques (JCGT) 3, 3, 7–34. http://jcgt.org/published/0003/03/02/.]が、我々は $\theta$ の最小値と最大値の算術平均に対する距離を単純に計算する。そして、以下のように透過率で重み付けされたPrefix Sumを計算する。

\begin{align} b_{\sum}(u, v) &:= \frac{\sum_{j=0}^u w_{j,v} \cdot b(j, v)}{\sum_{j=0}^u w_{j,v}} \\ w_{u,v} &:= \left[ -\frac{1}{\sigma_t} \cdot \exp(-\sigma_t \cdot t) \right]_{(u - \frac{1}{2})\Delta(v)}^{(u + \frac{1}{2})\Delta(v)} \end{align}

ある位置 $q \in \mathbb{R}^3$ で終わる光線の散乱を計算するため、適切な位置で事前フィルタされたシャドウマップ $b_{\sum}$ をサンプルし、1がフィルタされたハードシャドウになるように0から1の間でシャドウ強度を再構築し、可視性を得るために1から引いて、最大限可能な散乱で乗算する。

f(\omega_p, \omega_l) \cdot E_l \cdot \left[ -\frac{1}{\sigma_t} \cdot \exp(-\sigma_t \cdot t) \right]_{0}^{\|q - p\|^2}

この手続きはスクリーンでピクセルあたりに事前フィルタされたシャドウマップ $b_{\sum}$ における単一のルックアップのみを必要とする。故に、このテクニックの実行時間はシーンの複雑さに独立である。

5.2. フィルタリング付き修正（Rectification with Filtering）

[Klehm et al. 2014Klehm, O., Seidel, H.-P. and Eisemann, E. 2014. Filter-based real-time single scattering using rectified shadow maps. Journal of Computer Graphics Techniques (JCGT) 3, 3, 7–34. http://jcgt.org/published/0003/03/02/.]で提案される線形修正は、遠くのジオメトリが圧縮される一方で、カメラ近くのジオメトリに対してシャドウマップの大部分を割り当てる傾向にある。これはニアクリップ面を離したり、シャドウマップを分割することで緩和できるが、いずれのソリューションもまったく満足するものではない。加えて、非線形修正は依然としてepipoleが視野角に近いケースに対して実装されなければならない。

一方で、上記で述べられる非線形修正はラスタライゼーションハードシャドウで効率的に実装するためにシャドウマップのリサンプリングを必要とする。一般的なシャドウマップはリサンプリング中にフィルタできないので、大量のエイリアシングアーティファクトを引き起こす。直線的な輪郭は薄明光線《crepuscular ray》での｜明らかな《visible》縞模様を引き起こす階段状のアーティファクトを示す。これらの縞模様はカメラの移動や回転に関して安定ではないので、かなり目立ってしまう。

我々はMoment Shadow Mapを用いて、リサンプリング中のフィルタリングを達成した。非線形な方法で深度を歪ませるので、取得したモーメントを深度分布に変換し直す。小さなフィルタ領域で動作しているので、単純な分布を予期するのは適切である。ほとんどの場合、そのフィルタ領域は2つより多い表面をカバーしない。

節2.2.1では、 $z_0 = z_f$ を定める4つのモーメントから3つの深度値 $z_0, z_1, z_2$ を持つ深度分布を再構成する方法を示した。任意の選択を回避し、より効率的なソリューションを得るため、 $z_0$ を無限大に飛ばすとすると、 $w_0$ は0に近づき、深度値 $z_0$ を捨てることができる。残りの分布 $w_1 \delta_{z_1} + w_2 \delta_{z_2}$ はモーメント $b_0, b_1, b_2, b_3$ と未だに互換性がある。ほぼ定数の深度での2つ以下の表面の仮定の下、再構成が十分であること確信できる。そうでなければ、単一のシャドウマップサンプルより確実に良い、理に適った近似を提供する。

アルゴリズム2: 3つのモーメントから2つの深度値を持つ深度分布の構築
入力: $b_0 = 1$ と $b_2 - b_1^2 > 0$ を持つモーメントのベクトル $b \in \mathbb{R}^4$
出力: すべての $j \in \{0, 1, 2, 3\}$ に対する $\varepsilon_Z(\boldsymbol{z}^j) = b_j$ のような $\mathbb{R}$ 上の分布 $Z$
1. $q_2 := b_2 - b_1^2$ をセットする
2. $q_1 := b_1 \cdot b_2 - b_3$ をセットする
3. $q_0 := -b_1 \cdot q_1 - b_2 \cdot q_2$ をセットする
4. 解 $z1, z2 \in \mathbb{R}$ を得るために $q_2 \cdot z^2 + q_1 \cdot z + q_0 = 0$ を解く
5. $w_2 := \frac{b_1 - z_1}{z_2 - z_1}$ と $w_1 := 1 - w_2$ をセットする
6. $w_1 \cdot \delta_{z_1} + w_2 \cdot \delta_{z_2}$

述べられた分布はアルゴリズム2で構築される。これは非負の分散 $\sigma^2 := q_2 = b_2 - b_1^2$ の入力で失敗するが、この場合は $\delta_{b_1}$ を単純に戻すことで適切に扱われる。

分布を得たらば、節5.1.1に述べたように、その深度 $z_1, z_2$ を傾斜角 $\theta_1, \theta_2$ に変換して、範囲外の値をクランプして区間 $[-1, 1]$ に正規化する。そして、 $\Theta_m^* \circ \boldsymbol{b}$ を経由してモーメントのベクトルの両方の値を変換して、重み $w_1, w_2$ を用いて線形に組み合わせ、その結果は $b(u, v)$ に格納される。この時点で、我々は4つ以上のモーメントや他のフィルタ可能なシャドウマップを生成することもできる。

このスキームの使用は全体的に任意である。我々の実装はピクセルシェーダで修正されたフィルタ可能なシャドウマップ $b(u, v)$ を生成する。フィルタリングが有効である場合、ピクセルシェーダはMoment Shadow Mapからフィルタされたサンプルを取り、そうでなければ、通常のシャドウマップからテクセルを単に読み込む。Moment Shadow Mapからのサンプルは若干高価であるのみである(節5.5.2参照)。

5.3. 4つのモーメントを持つ事前フィルタされた単一散乱（Prefiltered Single Scattering with Four Moments）

事前フィルタされた単一散乱では、Convolution Shadow Mapの代わりに通常の量子化変換を持つMoment Shadow Mapを用いることができる。

ただし、単一散乱では一般に、サーフェスアクネを回避するためのシャドウ強度の過小推定を必要としない。我々のソリューションでは鋭い下界と上界の重み付けされた組み合わせを取る。上界はアルゴリズム1で下界に $w_0$ を加えることで得られる。故に、両界を計算するオーバーヘッドは小さい。

5.3.1. 適応的な過大推定（Adaptive Overestimation）

鋭い上界と下界を持つと、2つの間を補間する重み $\beta \in [0, 1]$ が必要になる。単一散乱は $\beta = 1$ で過小推定し、 $\beta = 0$ で過大推定するので、単純なアプローチではワーストケースの誤差が最小になるよう $\beta = \frac{1}{2}$ とするだろうが、その重みはピクセルあたりに任意にセットでき、より洗練された選択がアーティファクトを回避する。

$\beta$ が適切でないと、ライトリークが起こったり、フラグメントの深度値より大きい深度値がなくなったりする。

中間のケースで最適な近似品質を維持しつつ両方のアーティファクトを回避するため、 $\beta$ を視線の方向 $\omega_p$ に依存するようセットする。この重み $\beta$ は $\omega_p = -\omega_l$ のときに $0$ になり、 $\omega_p = \omega_l$ のときに $1$ になり、epipole近くでは緩やかに変化すると思われる。我々の実験では、以下の選択が滑らかでもっともらしい結果をもたらしつつ、上で述べたアーティファクトを確実に除去することを発見した。

\beta = \frac{1 + \omega_l^\mathsf{T} \cdot \omega_p}{2}

5.4. 6モーメントを持つ事前フィルタされた単一散乱（Prefiltered Single Scattering with Six Moments）

5.4.1. ルートの計算（Computation of Roots）

4x4線形システム $B(b) \cdot q = (z^0, \dots, z^3)^\mathsf{T}$ をコレスキー分解を用いて解くことは依然として上手く動作する。次に、三次方程式 $\sum_{j=0}^3 q_j \cdot z^j = 0$ を $z$ に対して解く必要がある。様々な反復的および閉形式の解を実験した後、 Blinn [2007Blinn, J. F. 2007. How to solve a cubic equation, part 5: Back to numerics. IEEE Computer Graphics and Applications 27, 3, 78–89. 10.1109/MCG.2007.60. https://courses.cs.washington.edu/courses/cse590b/13au/lecture_notes/solvecubic_p5.pdf.] により提案される閉形式の解の一種に決めた。

その論文に示されるアルゴリズムは、打ち消しを回避するための最大最小の大きさのルートの計算のための2つの異なる分岐を用いる。我々のアプリケーションでは、このオーバーヘッドが不必要であることが判明した。3つのルートすべてを計算するのに2の分岐の内の1つを用いると、アーティファクトフリーな結果を生み出す。他の閉形式の解は $|q_3| \ll 1$ でアーティファクトに悩まされ、反復的な解は計算オーバーヘッドが高かった。試みた解の中で、Blinnの研究に基づくものがもっとも高速である。

5.4.2. 境界の計算（Computation of Bounds）

最終段階では単一散乱の強さに比例する平均可視性を近似する。これはモーメント制約に従う重み $w_0, \dots, w_3$ の線形結合である。

\left( \begin{array}{c} 1 & 1 & 1 & 1 \\ z_0 & z_1 & z_2 & z_3 \\ z_0^2 & z_1^2 & z_2^2 & z_3^2 \\ z_0^3 & z_1^3 & z_2^3 & z_3^3 \end{array} \right) \cdot \left( \begin{array}{c} w_0 \\ w_1 \\ w_2 \\ w_3 \end{array} \right) = \left( \begin{array}{c} b_0 \\ b_1 \\ b_2 \\ b_3 \end{array} \right)

線形結合での重みは以下のようになる。

v_0 := \beta \text{ and } \forall l \in \{1,2,3\} : v_l := \begin{cases} 0 & \text{if } z_f > z_l \\ 1 & \text{if } z_f \le z_l \end{cases}

内積として書くと以下のようになる。

\left( \begin{array}{c} w_0 \\ w_1 \\ w_2 \\ w_3 \end{array} \right)^\mathsf{T} \cdot \left( \begin{array}{c} v_0 \\ v_1 \\ v_2 \\ v_3 \end{array} \right) = \left( \begin{array}{c} b_0 \\ b_1 \\ b_2 \\ b_3 \end{array} \right)^\mathsf{T} \cdot \underbrace{ \left( \begin{array}{c} 1 & z_0 & z_0^2 & z_0^3 \\ 1 & z_1 & z_1^2 & z_1^3 \\ 1 & z_2 & z_2^2 & z_2^3 \\ 1 & z_3 & z_3^2 & z_3^3 \end{array} \right)^{-1} \cdot \left( \begin{array}{c} v_0 \\ v_1 \\ v_2 \\ v_3 \end{array} \right) }_{=: u}

この最後の数式の行列はVandermonde行列であるので、ベクトル $u \in \mathbb{R}^4$ は $l \in \{0,1,2,3\}$ となる $z = z_l$ に対する $v_l$ を取る補間多項式 $\sum_{j=0}^3 u_j \cdot z^j$ の係数を保持する。差商《divided differences》を用いてその｜ニュートン形式《Newton form》を構築し、その後、多項式の標準基底《canonical basis》に変換し直す[Greenbaum and Chartier 2012, p.181 ff.Greenbaum, A. and Chartier, T. P. 2012. Numerical methods: Design, analysis, and computer implementation of algorithms.]。これは効率的に、我々の目的で十分ロバストに動作する。

5.4.3. 量子化とバイアス付け（Quantization and Biasing）

複雑な深度分布があることで、すべての計算を完璧に正確に行っても完全に再構成できない。故に、単一散乱に対してそれほど極端ではない強いバイアスの効果が期待され、低精度のMoment Shadow Mapを用いることは魅力的である。

再び、丸め誤差はMoment Shadow Mapに格納される前にモーメントに適用されるアフィン変換の方法によって減らされるべきである。節2.4にある通り、これを定めるために数値的な最適化を用いる。我々は一般的な変換や偶奇番目でモーメントを別々に処理する変換を実験した。見つかった最良の変換は後者のカテゴリに属す。これは以下で与えられ、格納されるデータのエントロピーをテクセルあたり12.5ビット増加させる。

\Theta_6^* \left( \begin{array}{c} b_1 \\ b_2 \\ b_3 \\ b_4 \\ b_5 \\ b_6 \end{array} \right) = \left( \begin{array}{c} \left( \begin{array}{c} 2.5 & -1.8749986 & 1.26583039 \\ -10 & 4.20757543 & -1.47644883 \\ 8 & -1.83257679 & 0.71061660 \end{array} \right)^\mathsf{T} \cdot \left( \begin{array}{c} b_1 \\ b_3 \\ b_5 \end{array} \right) + \left( \begin{array}{c} 0.5 \\ 0.5 \\ 0.5 \end{array} \right) \\ \left( \begin{array}{rrr} 4 & 9 & 0.57759806 \\ -4 & -24 & 4.61936648 \\ 0 & 16 & -3.07953907 \end{array} \right)^\mathsf{T} \cdot \left( \begin{array}{c} b_2 \\ b_4 \\ b_6 \end{array} \right) + \left( \begin{array}{r} 0 \\ 0 \\ 0.018888946 \end{array} \right) \end{array} \right)

最適なバイアスは4モーメントの場合と同じように決定される(節2.5を参照)。そのバイアスするモーメントのベクトルは以下となる。

b^* := (1, 0, 0.5566, 0, 0.489, 0, 0.47869382)^\mathsf{T}

我々のDirect3D 11の実装におけるモーメントのストレージは、10ビットの固定小数点数を格納する各3チャネルを持つ2つのテクスチャ(テクセルあたり64ビット、 $\alpha_b = 4 \cdot 10^{-3}$ )、16ビットの固定小数点数を格納する2チャネルと4チャネルのもの(テクセルあたり96ビット、 $\alpha_b = 8 \cdot 10^{-5}$ )、単精度floatを格納する3チャネルのもの2つ(テクセルあたり192ビット、 $\alpha_b = 5 \cdot 10^{-6}$ )のいずれかを用いる。

5.5. 結果と議論（Results and Discussion）

5.5.1. 質的評価（Qualitative Evaluation）

5.5.2. 実行時間（Run Time）

5.5.3. 結論（Conclusions）

6. おわりに（Conclusions）

謝辞（Acknowledgments）

A.

Footnotes

いわゆるアルファブレンディング。 $\alpha c_{src} + (1 - \alpha) c_{dest}$ ↩
$a+b$ の $a$ の方 ↩