このページではMathJaxによる数式表示を行っています。対応しているブラウザーであれば少し待てば数式が表示されます。
(画面左下にMathJaxのタイプセットステータスが表示されます。)
環境によっては表示が崩れているかもしれません。その際はブラウザーを変更してみてください。

関与媒質中の自由行程サンプリングと透過率の推定

みなさんこんにちは!shockerです。このページはレイトレ合宿5‽のアドベントカレンダー1週目の記事として制作しました。レイトレ合宿4!?の八丈島はやばい台風の予報のなか奇跡的に晴れましたが、2!!や3!!!は結構曇っててレイが飛びにくそうとかどうでもいいこと言ってましたね!今回はその曇りの原因である雲を始めとする関与媒質のレンダリングに必要な要素について紹介したいと思います。

ボリュームレンダリングで光の経路を構築するには、サーフェスレンダリングで用いた方向のサンプリングに加えて関与媒質中を光が媒質や物体表面と反応せずに進む距離、つまり自由行程のサンプリングが必要となります。また、与えられた距離に対する媒質の透過率の推定も必要となります。このページでは一様媒質・非一様媒質中それぞれの自由行程のサンプリング手法と透過率の推定手法について紹介します。

扱う問題

\begin{eqnarray*} L_i \Prta{\vx, \vomega} = T \, \Prta{\vx_S, \vx} L_o \Prta{\vx_S, \vomega} + \int_0^{\norma{\vx_S - \vx}} T \, \Prta{\vx', \vx} \Brk{ L_e^V \Prta{\vx', \vomega} + \sigma_s \Prta{\vx'} \int_{\mathcal{S}^2} f_p \Prta{\vx', \vomega', \vomega} L \Prta{\vx', \vomega'} d\vomega' } ds \label{eq:VLTE} \end{eqnarray*}

関与媒質をレンダリングする際には上記のボリュームレンダリング方程式を解くことになります。パストレーシングを始めとするモンテカルロ積分を用いる多くの手法では、光の輸送経路を逐次構築するために上式の積分に関してサンプリングを行う必要があります。内側の積分に関しては位相関数のサンプリング、つまり光の入射方向のサンプリングとなり、真空を仮定したサーフェスレンダリングの場合にも同様の概念がありました。しかしボリュームレンダリングの場合には、それに加えて外側の積分に関するサンプリングが概念的に増えています。これは現在の点 $ \vx $ から次の散乱点である $ \vx' $ もしくは $ \vx_S $ までの距離をサンプリングするという意味で、物理的な解釈としては自由行程(free path)のサンプリングに相当します。積分中の大括弧内や、$ L_o(\vx_S, \vomega) $ の項は再帰的に(ボリューム)レンダリング方程式を解かないと推定できないため、基本的な戦略としては透過率 $ T(\vx', \vx) $ に対応した重点的サンプリングを行うことになります。
物体表面からの輝度 $ L_o(\vx_s, \vomega) $ と媒質中からの輝度(積分中の大括弧内)を分けて推定すると、経路の分岐が爆発してしまうため、自由行程のサンプリング区間を積分区間 $ [0, \norma{\vx_S - \vx}) $ とするのではなく$ [0, \infty) $ として、自由行程が物体表面までの距離を超えた場合は、物体表面からの輝度を評価するという実装がよく使用されます。モンテカルロ推定関数を式で表すと次のようになります。

\begin{eqnarray} L_i \Prta{\vx, \vomega} \approx \begin{cases} \displaystyle \frac{T \, \Prta{\vx', \vx}}{p_d(s)} \Brk{ L_e^V \Prta{\vx', \vomega} + \sigma_s \Prta{\vx'} \int_{\mathcal{S}^2} f_p \Prta{\vx', \vomega', \vomega} L \Prta{\vx', \vomega'} d\vomega' } & s < \norma{\vx_S - \vx} \\ \displaystyle \frac{T \, \Prta{\vx_S, \vx}}{P_S(s \ge \norma{\vx_S - \vx})} L_o \Prta{\vx_S, \vomega} & s \ge \norma{\vx_S - \vx} \end{cases} \label{eq:estimator} \end{eqnarray}

ここで、$ \vx' = \vx - s \vomega $ です。2つの場合分けを確率的に選ぶようになるので確率密度の積分 + 確率は1になります。つまり、

\begin{equation*} \int_0^{\norma{\vx_S - \vx}} p_d(s) \, ds + P_S(s \ge \norma{\vx_S - \vx}) = 1 \end{equation*}

を満たす必要があります。

Next Event Estimation双方向パストレーシングにおけるサブパスの接続では、自由行程のサンプリングは行われず、与えられた2点間の透過率 $ T $ の計算(推定)を行う必要があります。(いずれのリンク先のページも、真空を前提としたアルゴリズムの説明となっているので透過率 $ T $ の項は登場しませんが、関与媒質を扱う場合には現れます。)

自由行程のサンプリング

確率密度関数 (PDF) と累積分布関数 (CDF)

前節では、自由行程のサンプリング区間を $ [0, \infty) $ として透過率 $ T $ に対応した重点的サンプリングを行うと書きました。これをそのままPDFが満たすべき条件式として書くと次のようになります。

\begin{equation*} \int_0^\infty p_d(s) \, ds = k \int_0^\infty T(s) \, ds = k \int_0^\infty \exp\Prtc{-\int_0^s \sigma_e(s') ds'} ds = 1 \end{equation*}

ここで、$ k $ は確率密度関数として正規化するための係数です。また、簡単のため $ T(s) = T(\vx, \vx') $ として書いています。係数 $ k $ やCDFを求めるためにPDFを積分する必要があります。$ \sigma_e(s') $ が定数であれば簡単そうですが、定数でない場合、つまり非一様媒質の場合にはかなり難しそうです。
ところで、透過率 $ T $ は $ s = 0 $ のとき $ 1 $ で $ s \to \infty $ のとき $ 0 $ です。一方求めたいCDFは $ s = 0 $ のとき $ 0 $ で $ s \to \infty $ のとき $ 1 $ となります。この関係より、透過率 $ T $ が既にほとんどCDFのようなものを表していると気がつけば、PDF, CDFの導出はかなり簡単になります。したがって、CDFは次の式で表されます。

\begin{equation} P_d(s) = 1 - T(s) = 1 - \exp\Prtc{-\int_0^s \sigma_e(s') ds'} \label{eq:CDF} \end{equation}

PDFはCDFを微分することで求まるので次のようになります。

\begin{equation} p_d(s) = \frac{d}{ds} P_d(s) = \sigma_e(s) \exp\Prtc{-\int_0^s \sigma_e(s') ds'} = \sigma_e(s) T(s) \label{eq:PDF} \end{equation}

これで目的のPDFとCDFが求まりました。注意点としては、このPDFは$ \sigma_e(s) $ が定数ではない場合、透過率だけに厳密に比例していないことが挙げられます。しかし、式 $ \eqref{eq:estimator} $ の媒質中の場合に代入してみるとわかりますが、透過率以外の項に含まれる媒質に関する係数も同時に考えれば、むしろこのPDFのほうが好ましくなります。
なお、自由行程が物体表面までの距離を超える確率( = 物体表面からの輝度を評価する確率)はCDFを使って求めることができて、次の式のようになります。

\begin{equation} P_S(s \ge \norma{\vx_S - \vx}) = 1 - P_d(\norma{\vx_S - \vx}) = \exp\Prtc{-\int_0^{\norma{\vx_S - \vx}} \sigma_e(s') ds'} = T(\norma{\vx_S - \vx}) \label{eq:prob_surface} \end{equation}

以上により、モンテカルロ推定関数 $ \eqref{eq:estimator} $ に必要なPDF(と確率)の式を得ることができました。が、結局これらの式の解析的な値を知ることは一般的なケースにおいて不可能です。しかし幸いなことに推定関数の分子にある透過率とPDF(or 確率)が解析的にわからない積分部をキャンセルし合います。つまり、被積分関数と確率密度関数それぞれの値を明示的に知ることは不可能ですが、それらから計算されるモンテカルロ推定関数(のウェイト)の値は知ることができます。
PDF $ \eqref{eq:PDF} $ と確率 $ \eqref{eq:prob_surface} $ を用いた場合のモンテカルロ推定関数は次のようになります。

\begin{eqnarray*} L_i \Prta{\vx, \vomega} \approx \begin{cases} \displaystyle \frac{1}{\sigma_e(\vx')} \Brk{ L_e^V \Prta{\vx', \vomega} + \sigma_s \Prta{\vx'} \int_{\mathcal{S}^2} f_p \Prta{\vx', \vomega', \vomega} L \Prta{\vx', \vomega'} d\vomega' } & s < \norma{\vx_S - \vx} \\ \displaystyle \vphantom{\frac{T}{P_S}} L_o \Prta{\vx_S, \vomega} & s \ge \norma{\vx_S - \vx} \end{cases} \end{eqnarray*}

一様媒質中のサンプリング

一様媒質における自由行程のサンプリングは非常に簡単で、逆関数法を使うことができます。
CDFは 式 $ \eqref{eq:CDF} $ より次のようになります。

\begin{equation*} P_d(s) = 1 - \exp\Prta{-\sigma_e s} \end{equation*}

逆関数法を使うと、自由行程は次の式でサンプリングできます。

\begin{equation*} s = -\frac{\ln(1 - u)}{\sigma_e} \end{equation*}

ここで、$ u $ は $ [0, 1) $ 中に分布する一様乱数です。

非一様媒質中のサンプリング

非一様媒質の場合は、以下に示すように光学的深度に関する式までは変形できますが、それ以上の変形が難しく逆関数法を使うことができません。

\begin{equation} \int_0^s \sigma_e(s') ds' = -\ln(1 - u) \label{eq:sampling_heterogeneous} \end{equation}

そこで以下に示すような手法が自由行程のサンプリングのために使われます。

レイマーチング (Ray Marching)

図1. レイマーチングによる自由行程サンプリング
まず光学的深度をサンプリングし、その光学的深度 $ \tau_{sample} $ に到達するまで小区間ごとに積分していく。この図では0次近似だが、台形則(1次近似)を使えば連続的な自由行程サンプリングも行える。

関与媒質中のレイの位置を小さな区間ずつ進めながら式 $ \eqref{eq:sampling_heterogeneous} $ を満たす距離を探す手法をレイマーチング(ray marching)と呼びます(図1)。基本的には各小区間内は消散係数が定数であるという近似の下で式 $ \eqref{eq:sampling_heterogeneous} $ の積分を計算するため、実現される自由行程の分布は厳密に式 $ \eqref{eq:PDF} $ には沿わず、小区間の幅に応じてバイアスのかかったサンプリングになります。
以下にレイマーチングによる自由行程サンプリングの疑似コードを示します。

function $ \mathrm{rayMarching} $(レイの原点 $ \vx $, 入射方向 $ \vomega $, 初期距離 $ s_0 $, 媒質境界までの距離 $ s_e $) {
	光学的深度 $ \tau_{sample} $ をサンプル $ \tau_{sample} = -\ln(1 - u) $
	自由行程を初期化 $ s \LAR s_0 $
	現在の光学的深度 $ \tau $ を初期化 $ \tau = 0 $
	while $ \tau < \tau_{sample} $ {
		現在の位置 $ \vx' = \vx - s \cdot \vomega $
		現在の位置の消散係数を得る $ \sigma_e(\vx') $
		光学的深度を更新 $ \tau \LAR \tau + \sigma_e(\vx') \cdot \Delta s $
		位置を小区間 $ \Delta s $ 進める $ s \LAR s + \Delta s $
		if 現在の位置が関与媒質の外 $ s > s_e $
			break
	}
	return $ s $
}

デルタトラッキング (Delta Tracking)

図2. デルタトラッキングによる自由行程サンプリング
一様化した消散係数 $ \bar{\sigma}_e $ を使って自由行程をサンプリングし、仮の散乱点を生成、それが本来の媒質による散乱かどうかを確率的に判定し、偽の散乱であれば自由行程を延長するという処理を繰り返す。オレンジ色の短い線がそれぞれの仮の散乱点において使用された一様乱数を表している。

偽の粒子で関与媒質を埋めることによって消散係数を一様化(homogenize)し、一様媒質と同様な自由行程のサンプリングを可能とする手法をデルタトラッキング(delta tracking, Woodcock tracking、別の分野ではさらにpseudo scattering, null-collision algorithm)と呼びます。ただし、サンプリングされた自由行程が適当かどうかを散乱点の媒質特性に応じて確率的に判定することにより、非一様媒質としての正しいサンプリングを実現します(図2)。
具体的な手順としては、まず一様化された消散係数 $ \bar{\sigma}_e $ を用いて一様媒質と同様の解析的なサンプリングを行い「仮の」自由行程を得ます。その後、仮の散乱点における本来の消散係数 $ \sigma_e(s) $ の割合 $ \sigma_e(s) \, / \, \bar{\sigma}_e $ を評価します。この割合に対応した確率で実際の粒子か偽の粒子いずれかとの散乱が発生します。偽の粒子の散乱アルベドは $ 1 $、位相関数は $ \delta(\vomega) $ として設定することで、偽の粒子との散乱は何の散乱も起こしていないことに等しくなります。偽の粒子との散乱と判定された場合には、再度 $ \bar{\sigma}_e $ を用いて自由行程を延長し、同様の処理を繰り返します。この一連の処理によって式 $ \eqref{eq:PDF} $ に厳密に沿ったバイアスの無い自由行程のサンプリングが可能になります。
なお、本来の消散係数の割合をそのまま有効な確率として使うために、一様化された消散係数 $ \bar{\sigma}_e $ は自由行程サンプリングを行う直線上において常に本来の消散係数 $ \sigma_e(s) $ 以上の値を持つ必要があります。

\begin{equation*} \bar{\sigma}_e \ge \sigma_e(s) \end{equation*}

この理由により一様化された消散係数をmajorant extinction coefficient、短くはmajorantと呼びます。majorantはいくらでも大きな値に設定することができますが、大きなmajorantは短い自由行程の延長距離を意味し、計算コストの増加につながります。最低限のmajorantにおいてデルタトラッキングは最も効率が良くなります。
以下にデルタトラッキングによる自由行程サンプリングの疑似コードを示します。

function $ \mathrm{deltaTracking} $(レイの原点 $ \vx $, 入射方向 $ \vomega $, 初期距離 $ s_0 $, 媒質境界までの距離 $ s_e $) {
	$ [0, 1) $ の一様乱数 $ u_0 $ を用いて自由行程をサンプリング $ \displaystyle s \LAR s_0 - \frac{\ln(1 - u_0)}{\bar{\sigma}_e} $
	現在の位置 $ \vx' = \vx - s \cdot \vomega $
	if 現在の位置が関与媒質の外 $ s > s_e $
		return $ s $
	while $ [0, 1) $ の一様乱数 $ u_1 $ に関して $ \displaystyle u_1 \ge \frac{\sigma_e(\vx')}{\bar{\sigma}_e} $ {
		$ [0, 1) $ の一様乱数 $ u_0 $ を用いて自由行程を延長 $ \displaystyle s \LAR s - \frac{\ln(1 - u_0)}{\bar{\sigma}_e} $
		現在の位置 $ \vx' = \vx - s \cdot \vomega $
		if 現在の位置が関与媒質の外 $ s > s_e $
			break
	}
	return $ s $
}

透過率の推定

一様媒質中の推定

一様媒質における透過率の推定は単純に次の式を計算するだけです。確率的な要素も絡まないので推定というより、もはやただの計算です。

\begin{equation*} T(\vx', \vx) = \exp\Prta{-\sigma_e \norma{\vx' - \vx}} \end{equation*}

非一様媒質中の推定

非一様媒質における透過率の式は次の形で表されます。

\begin{equation} T(\vx', \vx) = \exp\Prtc{-\int_0^{\norma{\vx' - \vx}} \sigma_e(s') ds'} \label{eq:transmittance} \end{equation}

この式から解析的な値を得ることはできないので、何かしらの数値計算手法で透過率を推定する必要があります。

指数関数 $ \exp $ によるバイアス

上の式 $ \eqref{eq:transmittance} $ を見て、「素直にモンテカルロ法とか使って推定してしまえるんじゃない?」と思うかもしれませんが、ここには罠があります。光学的深度(積分部)をモンテカルロ法でバイアス無く推定できたとしても、透過率はその値を指数関数 $ \exp $ にかけたものであり、その期待値は真値からずれてしまいます。
光学的深度の推定値を確率変数 $ X $、そこから計算される透過率の推定値を確率変数 $ Y = \exp(-X) $ とおきます。期待値のずれは以下のように表すことができます。

\begin{equation*} \exp(-E[X]) \neq E[Y] = E[\exp(-X)] \end{equation*}

光学的深度の推定がバイアス無く行えている場合は左辺が真値となります。ただし、それを達成するためには光学的深度に関して無限のサンプル数が必要です。そんなことをしていては、パストレーシングなどのアルゴリズムの1サンプルに無限に時間がかかってしまうので、実際には有限のサンプル数で推定した光学的深度から透過率を計算することになりますが、達成される推定は右辺のようになってしまい、期待値が真値からはずれてしまいます。

レイマーチング (Ray Marching)

透過率の推定にもレイマーチングを使うことができます。自由行程のサンプリングと同様に近似が入るために推定結果にはバイアスがかかります。各小区間中でサンプル点をランダムに選べば、光学的深度の推定はバイアス無く行えます[Pauly1999]。しかし前述の理由により透過率の推定としては結局バイアスがかかってしまいます。

function $ \mathrm{rayMarchingEstimator} $(レイの原点 $ \vx $, 入射方向 $ \vomega $, 一端の距離 $ s_0 $, もう一端の距離 $ s_1 $) {
	光学的深度 $ \tau $ を初期化 $ \tau = 0 $
	小区間の幅を計算 $ \displaystyle \Delta s \LAR \frac{s_1 - s_0}{N} $
	for $ i \RAR 0 $ to $ N - 1 $ {
		現在の距離 $ s \LAR s_0 + i \cdot \Delta s $
		現在の位置 $ \vx' = \vx - s \cdot \vomega $
		現在の位置の消散係数を得る $ \sigma_e(\vx') $
		光学的深度を累積 $ \tau \LAR \tau + \sigma_e(\vx') \cdot \Delta s $
	}
	return $ \exp(-\tau) $
}

デルタトラッキング (Delta Tracking)

デルタトラッキングは透過率の推定に使うことも可能で、バイアスも無くせます。デルタトラッキングにより生成された自由行程が、透過率を計算する区間より長かった場合は $ 1 $、短かった場合は $ 0 $ という二値の確率変数になります。

function $ \mathrm{deltaTrackingEstimator} $(レイの原点 $ \vx $, 入射方向 $ \vomega $, 一端の距離 $ s_0 $, もう一端の距離 $ s_1 $) {
	デルタトラッキングにより自由行程を得る $ s = \mathrm{deltaTracking} $($ \vx $, $ \vomega $, $ s_0 $, $ s_1 $)
	return $ s > s_1 $
}

レシオトラッキング (Ratio Tracking)

デルタトラッキングでは仮の散乱点における本来の消散係数の割合、比に応じた確率的な自由行程の延長を行ない、透過率の推定では最終的に二値の値を返します。確率的な処理の打ち切りの代わりに、結合確率を保持・更新し連続的な値を返すように変更を加えた手法をレシオトラッキング(ratio tracking)[Novák2014]と呼びます。レシオトラッキングでは、結果の分散は減少しますが端点に到達するまで必ず処理を継続するため計算負荷は高くなります。

function $ \mathrm{ratioTrackingEstimator} $(レイの原点 $ \vx $, 入射方向 $ \vomega $, 一端の距離 $ s_0 $, もう一端の距離 $ s_1 $) {
	ウェイトを初期化 $ T = 1 $
	$ [0, 1) $ の一様乱数 $ u_0 $ を用いて自由行程をサンプリング $ \displaystyle s \LAR s_0 - \frac{\ln(1 - u_0)}{\bar{\sigma}_e} $
	現在の位置 $ \vx' = \vx - s \cdot \vomega $
	if 端点に到達 $ s > s_1 $
		return $ T $
	while true {
		ウェイトを更新 $ \displaystyle T \LAR T \times \Prt{1 - \frac{\sigma_e(\vx')}{\bar{\sigma}_e}} $
		$ [0, 1) $ の一様乱数 $ u_0 $ を用いて自由行程を延長 $ \displaystyle s \LAR s - \frac{\ln(1 - u_0)}{\bar{\sigma}_e} $
		現在の位置 $ \vx' = \vx - s \cdot \vomega $
		if 端点に到達 $ s > s_1 $
			return $ T $
	}
}

さいごに

このページでは、レイマーチングやデルタトラッキング、レシオトラッキングといった自由行程サンプリングや透過率の推定のための手法を紹介しました。これらは簡単な手法ですが、発展的な手法[Novák2014, Szirmay2017, Kutz2017]もこれらを基礎として開発されています。

参考文献

  1. [Kutz2017] Peter Kutz, Ralf Habel, Yining Karl Li, Jan Novák - "Spectral and Decomposition Tracking for Rendering Heterogeneous Volumes", 2017
  2. [Novák2014] Jan Novák, Andrew Selle, Wojciech Jarosz - "Residual Ratio Tracking for Estimating Attenuation in Participating Media", 2014
  3. [Szirmay2017] László Szirmay-Kalos, Iliyan Georgiev, Milán Magdics, Balázs Molnár, Dávid Légrády - "Unbiased Light Transport Estimators for Inhomogeneous Participating Media", 2017
  4. [Pauly1999] Mark Pauly - "Robust Monte Carlo Methods for Photorealistic Rendering of Volumetric Effects", 1999

Today : 00000071 Total : 00061010
Copyright © 2017 shocker.