こんにちは。Ray-Tracingはメモリ消費量、計算量、精度の観点で様々な手法が開発されてきました。その中でも、Radianceを事前計算しておき、球面調和関数の展開係数として保存しておくことでレンダリング計算量とメモリ消費量を抑えることができる手法を紹介する。
#目次
1.はじめに
2.アルゴリズム
3.半球面調和関数
[4.incoming radianceとBRDFの半球面調和関数展開](#4-incoming radianceとBRDFの半球面調和関数展開)
5.半球面調和関数の勾配
6.半球面調和関数の補間（移動まで）
7.次回に向けて
8.参考文献
#1. はじめに
Ray-Tracingは関節光の効果を反映させることができるためグラフィックスのリアリティが格段に良くなる。しかし、計算量が大変多いため、ゲームなどのリアルタイムレンダリングでは、実装することは困難である。そのため、精度の低下を抑えつつ、メモリ消費量や計算量を抑える手法が数多く考案されてきた。今回紹介する手法は、Hemispherical Harmonic Radiance Interpolationと呼ばれる手法である。元にした論文は、[GKPB04]であり、球面調和関数変換の理論と実装で用いた定式化を利用する。
上の図が、今回紹介する手法の概念図である。事前計算の段階で黒線で書かれた半球に入射するRadianceを計算しておき、半球面調和関数の展開係数（HSH係数）として保存（cache）しておきます。レンダリング時、レイがヒットした場所で赤半球に入射する輝度は黒半球に入射する輝度として保存しておいたデータを内挿補完することで計算される。
#2. アルゴリズム
//事前計算
for (cacheデータを生成する可能性のある面)
if(cacheデータを生成するのに適した面)
その位置の座標と法線情報を取得
半球面に入射するRadianceを計算しHSH係数を保存：➀
HSH係数の勾配を計算して保存：➁
BRDFのHSH係数を保存：➂
end if
end for
//Radiance Cachingを用いたレンダリング
for (レイがSurfaceにヒットするたび)
近傍のRadiance Cache RecordからHSH係数を取得する
利用するHSH係数に掛ける重みを計算
HSH係数の勾配から、変化量を計算する
HSH係数を回転させて、座標系を揃える
内挿を行う
内挿結果とBRDFより、Outgoing Radianceを計算する
end do
#3.半球面調和関数
Ray-Tracingを行う時、必要な情報はある面上の点の上半球から入射するRadianceや上半球のBRDFだけであることがほとんどである。そのため、サンプリングを行う時に無駄な計算が増えたり、上半球と下半球の間でノイズが生じたりするので、効率的に上半球のデータを取得することを目的として、半球面調和関数が考案された。半球面調和関数と球面調和関数の関係は球面調和関数変換の理論と実装で記している定式化を用いると、簡単に書ける。
注意：$\widetilde{A}$のようなチルダ付き記号は半球面調和関数用の記号であることを示す。今後、式番号は球面調和関数変換の理論と実装で使われたものと区別するために、H付きの番号を利用する。
\widetilde{P}^m_n (\cos \theta) = {P}^m_n (2\cos \theta -1) \tag{H1}
H_{nm}(\theta,\phi)\equiv
\begin{cases}
\sqrt{2}\tilde{X}_{n|m|}(\theta) ~\cos (m\phi)~~~~(-n<m<0)\\
\tilde{X}_{n|m|}(\theta)~~~~~~~~~~~~~~~~~~~~~~~~~(m=0)\\
\sqrt{2}\tilde{X}_{nm}(\theta) ~\sin(m\phi)~~~~(-n<m<0) \tag{H2}
\end{cases}
\begin{eqnarray*}
\tilde{K}_{nm} &=& \biggl[ \frac{2n+1}{2\pi}\frac{(n-|m|)!}{(n+|m|)!} \biggr]^{\frac{1}{2}}\\
&=& \sqrt{2} K_{nm}
\end{eqnarray*}
\begin{eqnarray*}
\tilde{X}_{nm} (\theta) &=& \tilde{K}_{nm} \tilde{P}^{|m|}_n(\cos\theta)\tag{H3}\\
&=& \sqrt{2} X_{nm} (\theta)
\end{eqnarray*}
可視化すると以下のような形状として表現される。球面調和関数を上半球に寄せたような見た目となる。
#4.Incoming RadianceとBRDFの半球面調和関数展開
事前計算として、HSH係数のデータセットを作成しておく必要がある。半球面調和関数展開と球面調和関数展開のやり方はほとんど同じである。以下に、その過程を記す。
➀上半球に一様なサンプル$(\theta_i,\phi_i)$を生成する。球面上での一様乱数生成を参考にしてもらうとできる。上半球側だけに乱数を生成したいので、下半球側に生成されたサンプルはZ=0の平面を基準にして面対称な移動をしてやれば良い。
➁Incoming RadianceとBRDFをHSH係数に変換する。
\hat{\lambda}_{nm} = \frac{2\pi}{MN} \sum^{MN}_{i=1} L_i(\theta_i,\phi_i) H_{nm}(\theta_i,\phi_i) \tag{H4}
\hat{f}_{nm}(\theta_o,\phi_o) = \frac{2\pi}{MN} \sum^{MN}_{i=1} f_i(\theta_i,\phi_i,\theta_o,\phi_o) H_{nm}(\theta_i,\phi_i) \tag{H5}
#5.半球面調和関数の勾配
HSH係数の内挿補間を行うためには、以下に記すHSH係数の勾配を求める必要がある。この時参照にしている座標系はレイのヒット面に垂直なZ座標を持つローカルな座標系である。
\nabla \lambda_{nm} = \bigg[\frac{\partial \lambda_{nm}}{\partial x} ,\frac{\partial \lambda_{nm}}{\partial y},0\bigg] \tag{H6}
H6を具体的に計算する方法は数値的な方法と解析的な方法の2種類ある。
#####数値的手法
数値的な手法では、以下の図のように、対象とする点pの位置をdxだけずらすことで、$\partial \lambda_{nm} / \partial x = (\lambda'_{nm}-\lambda)/\Delta x$を計算する方法だ。点$p$を$p'$に移動することで、変わる量と変わらない量があることに注意する。図では赤色の量は変わらない。
記号の説明
$\Delta A_k$：pから$(\theta_k,\phi_k)$の方向にレイを飛ばして、初めに衝突するオブジェクトの微小面積
$n_k$：衝突面の法線
$\xi_k$：$n_k$とレイの間の角
$q_k$：衝突点の座標
$r_k = |p - q_k|$、 $r'_k = |p' - q_k|$：$p$や$p'$と$q_k$の距離
$\Omega_k = \Delta A_k \frac{\cos \xi_k}{r^2_k} = \frac{2\pi}{MN}$
\Omega'_k = \Delta A_k \frac{\cos \xi'_k}{r'^2_k} = \frac{2\pi}{NM} \frac{r^2_k}{r'^2_k}\frac{\cos \xi'_k}{\cos \xi_k} \tag{H7}
これらを用いると、$\lambda'_{nm}$は以下のように書ける。
\begin{eqnarray*}
\hat{\lambda'}_{nm} &=& \sum^{MN}_{k=1} \Omega'_k L_i(\theta_k,\phi_k) H_{nm}(\theta'_k,\phi'_k) \\
&=& \frac{2\pi}{MN}\sum^{MN}_{k=1} \frac{r^2_k}{r'^2_k}\frac{\cos \xi'_k}{\cos \xi_k} L_i(\theta_k,\phi_k) H_{nm}(\theta'_k,\phi'_k) \tag{H8}
\end{eqnarray*}
ここで、$L_i(\theta_k,\phi_k)$は移動前後で不変としている点に注意する必要があります。BRDFがlambertでない限りは等しいとは言えませんが、dxが十分小さいと考え、$L_i(\theta_k,\phi_k)=L_i(\theta'_k,\phi'_k)$としています。以上のことと、y方向にも同様の議論を行うと、H6が計算できます。
#####解析的手法
数値的な手法と同様に、$L_i(\theta_k,\phi_k)=L_i(\theta'_k,\phi'_k)$という仮定を置きます。そしてこの手法では、実際にH4をxで微分する所から始めます。
\begin{eqnarray*}
\frac{\partial \hat{\lambda}_{nm}}{\partial x} &=& \sum^{MN}_{k=1} \frac{\partial}{\partial x} \bigg(\Omega_k L_i(\theta_k,\phi_k) H_{nm}(\theta_k,\phi_k)\bigg)\\
&=& \sum^{MN}_{k=1} L_i(\theta_k,\phi_k) \bigg(\frac{\partial \Omega_k} {\partial x}H_{nm}(\theta_k,\phi_k) + \Omega_k\frac{\partial H_{nm}(\theta_k,\phi_k)}{\partial x} \bigg) \tag{H9}
\end{eqnarray*}
以下に、H9を計算するための式群を記します。一部の証明は別の記事にて紹介します。また、利用する文字は数値的手法で掲載している図と同じものを用いています。
H9括弧内第一項で利用
\frac{\partial \Omega_k}{\partial x} = \frac{2\pi}{N}\frac{r_k n_x + 3q_x\cos\xi_k}{r^2_k \cos\xi_k} \tag{H10}
H9括弧内第二項で利用
\frac{\partial H_{nm}(\theta_k,\phi_k)}{\partial x} = \frac{\partial \theta_k }{\partial x}\frac{\partial H_{nm}(\theta_k,\phi_k)}{\partial \theta_k} + \frac{\partial \phi_k }{\partial x}\frac{\partial H_{nm}(\theta_k,\phi_k)}{\partial \phi_k} \tag{H11}
独立変数の変換
\frac{\partial \theta_k }{\partial x} = -\cos \theta_k \cos \phi_k/r_k \tag{H12}
\frac{\partial \phi_k }{\partial x} = \sin \phi_k/(r_k \sin \theta_k) \tag{H13}
\frac{\partial \theta_k }{\partial y} = -\cos \theta_k \sin \phi_k/r_k \tag{H14}
\frac{\partial \phi_k }{\partial y} = -\cos \phi_k/(r_k \sin \theta_k) \tag{H15}
半球面調和関数の微分
\frac{\partial H_{nm}(\theta,\phi)}{\partial \theta}\equiv
\begin{cases}
-2\sqrt{2}\tilde{K}_{nm}~\cos (m\phi)\sin\theta \frac{dP_{nm}}{dx}(2\cos\theta -1) ~~~~(m>0)\\
-2\sqrt{2}\tilde{K}_{nm}~\sin (m\phi)\sin\theta \frac{dP_{n|m|}}{dx}(2\cos\theta -1) ~~~~(m<0)\\
-2 \tilde{K}_{n0} \sin{\theta} \frac{dP_{n0}}{dx}(2\cos\theta -1)
~~~~~~~~~~~~~~~~~~~~~~~~~~~~(m=0)\tag{H16}
\end{cases}
\frac{\partial H_{nm}(\theta,\phi)}{\partial \phi}\equiv
\begin{cases}
-m H_{n,-m}(\theta,\phi) ~~~~(m\neq0)\\
0
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~(m=0)\tag{H17}
\end{cases}
ルジャンドル陪関数の微分
\frac{\partial P_{nm}}{\partial x}(x) \equiv
\begin{cases}
\frac{1}{x^2-1}\Big( xnP_{nm}(x) - (m+n)P_{n-1,m}(x)\Big)~~~~~(m<n)\\
-(-1)^m x (2m-1)!!(1-x^2)^{\frac{m}{2}-1} ~~~~~~~~~~~~~~~~(m=n)\tag{H18}
\end{cases}
H1～H3、H10～H18と球面調和関数変換の理論と実装で用いた式を利用すれば、H9が計算でき、同様の計算をすれば、H6が解析的に求められる。
ここまでの過程で求めた$\nabla \lambda_{nm}$をHSH係数として保存しておけば良い。
#6.半球面調和関数の補間（移動まで）
ここまでは事前にHSH係数を計算しておき、データを保存するまでの過程を記述した。ここからは、レンダリング計算時に毎フレーム行う計算を説明する。大きく分けると、３ステップで計算できる。
➀各HSH係数の重みを計算する。
➁勾配を用いて補間する。
➂座標軸を合わせるためにHSH係数に回転を施す。
今回の記事は移動➁までを説明するが、次回の記事ではその続きを説明する。
#####HSH係数の重み
対象とする点$p$のRadianceを取得するために、i番目のHSH係数の重みは以下のように計算できる。
w_i(p) = \frac{1}{\frac{|p-p_i|}{R_i} + \sqrt{1-n \cdot n_i}} \tag{H19}
$R_i$はharmonic mean distance。$p_i$がどれくらい遮蔽されているかの指標。場合によっては事前に計算しておくこともできる。
R_i = \frac{MN}{\sum^{MN}_{k=1}\frac{1}{r_k} } \tag{H20}
#####HSH係数の配列順序
HSH係数をベクトルとして一つのまとまりで考える方が便利であるため、以下のような順序で並べることとする。ただし、球面調和関数変換の理論と実装で説明した順番とは異なることに注意すべきである。
\Lambda = \Big[ \lambda_{0,0} ~,~ \lambda_{1,-1} ~,~\lambda_{1,0} ~,~ \lambda_{1,1} ~,~ \lambda_{2,-2}~,~\cdots \Big] \tag{H21}
\frac{\partial \Lambda}{\partial x} = \Big[ \frac{\partial \lambda_{0,0}}{\partial x} ~,~ \frac{\partial \lambda_{1,-1}}{\partial x} ~,~\frac{\partial \lambda_{1,0}}{\partial x} ~,~ \frac{\partial \lambda_{1,1}}{\partial x} ~,~ \frac{\partial \lambda_{2,-2}}{\partial x}~,~\cdots \Big] \tag{H22}
y方向の微分も同様に定義できる。
#####移動補間（translational gradient Interpolation）
$p-p_i = (d_x,d_y,d_z)$とする。ただし、鉛直(z)方向の変化は小さいと見なして、$d_z=0$とする。$p_i$でのHSH係数の$p$への移動は
\Lambda'_i = \Lambda_i + d_x \frac{\partial \Lambda_i}{\partial x} + d_y \frac{\partial \Lambda_i}{\partial y} \tag{H23}
#7.次回に向けて
今回の記事では、Hemispherical Harmonic Interpolationをするために、Radiance Cachingを生成することが、話の大部分を占めていた。そして、Interpolationの移動の話まで書けたが、回転以降の話は次回の記事に書くとします。
#8.参考文献
[GKPB04]
Pascal Gautron, Jaroslav Kˇrivánek, Sumanta N. Pattanaik, and Kadi Bouatouch. A novel hemispherical basis for accurate and efficient rendering. In Eurographics Symposium on Rendering (to appear), June 2004.
[KGPB05]
Jaroslav Kˇrivánek, Pascal Gautron,Sumanta N. Pattanaik, and Kadi Bouatouch 2005. Radiance caching for efficient global illumination computation. Transactions on Visualization and Computer Graphics