Introduction to volume rendering @yumcyawiz レイトレ合宿8 セミナー

Volume Renderingとは? ボリューム(霧, 水, 肌)などをレンダリングすること 物体表面だけではなく、3次元空間中で光の吸収・散乱が起こる ボリュームレンダリング方程式(VRE) を解くことでレンダリングが出来る 煙 表面化散乱

Volume Renderingの基礎概念 ボリュームレンダリングは以下の3要素から構成される. 発光 吸収 散乱 吸収の強さは 吸収係数 で表される. 散乱の強さは 散乱係数 で表される. ある点 から の間でどれだけ光が減衰するかは 透過率 で表される. μ ​ (x, λ) a μ ​ (x, λ) s x ​ 0 x ​ 1 T(x ​ , x ​ ) 0 1

放射輸送方程式 媒質中のある点 における放射輝度の変化は以下の 放射輸送方程式 によって記述される. : 点 において方向 から来る放射輝度 : 点 における吸収係数 : 点 における散乱係数 : 点 において方向 に発光する放射輝度 : 点 において方向 へ散乱される外から来た放射輝度 x L(x, ) ω x ω μ ​ (x) a x μ ​ (x) s x L ​ (x, ) e ω x ω L ​ (x, ) s ω x ω ( ⋅ ω ∇)L(x, ) = ω −μ ​ (x)L(x, ) − a ω μ ​ (x)L(x, ) + s ω μ ​ (x)L ​ (x, ) + a e ω μ ​ (x)L ​ (x, ) s s ω

Volume Rendering Equation(VRE) 放射輸送方程式を距離 で積分することで, 以下のVREが得られる. 発光 + 散乱光 + 媒質終端から入射する光 [0, t] ​ ​ L(x, ) ω = ​ T(x, x − s )(μ ​ (x − s )L ​ (x − s , ) + μ ​ (x − s )L ​ (x − s , ))ds ∫ 0 t ω a ω e ω ω s ω s ω ω + T(x, x − t )L(x − t , ) ω ω ω

VREのモンテカルロ積分 積分変数が距離 なので, をランダムにサンプリング(自由行程サンプリング)すると : 距離 の確率密度関数 : 指示関数 : 距離 が 以上の値を取る確率 ​ ​ L(x, ) ω = ​ T(x, x − s )(μ ​ (x − s )L ​ (x − s , ) + μ ​ (x − s )L ​ (x − s , ))ds ∫ 0 t ω a ω e ω ω s ω s ω ω + T(x, x − t )L(x − t , ) ω ω ω s s p ​ (s) S S 1 ​ A P(S ≥ t) S t ​ ​ Term1(s) Term2 (x, ) L ω = T(x, x − s )(μ ​ (x − s )L ​ (x − s , ) + μ ​ (x − s )L ​ (x − s , )) ω a ω e ω ω s ω s ω ω = T(x, x − t )L(x − t , ) ω ω ω = ​ 1 ​ + ​ 1 ​ p ​ (s) S Term1(s) {s

自由行程サンプリング Question 距離 のサンプリングにどのような確率分布を用いるべきか? 透過率 に注目し, それに比例するように重点的サンプリングを行いたい. ​ ​ L(x, ) ω = ​ T(x, x − s )(μ ​ (x − s )L ​ (x − s , ) + μ ​ (x − s )L ​ (x − s , ))ds ∫ 0 t ω a ω e ω ω s ω s ω ω + T(x, x − t )L(x − t , ) ω ω ω s → T(x, x − s ) ω

自由行程サンプリング 簡単のために, まずは吸収/散乱係数が に依存しない 場合(一様媒質)で考える. 透過率の定義は以下の通り. , は今, 定数なので 一様媒質の場合, 透過率は指数関数になる x μ (x) t T(x, x − s ) ω = μ ​ (x) + μ ​ (x) a s = e− ​ μ ​ (x−s )ds ∫0 s t ′ ω ′ μa μs T(x, x − s ) ω = e− ​ μ ​ (x−s )ds ∫0 s t ′ ω ′ = e− ​ μ ​ ds ∫0 s t ′ = e−μ ​ s t

自由行程サンプリング 透過率は指数関数なので, 距離 の確率分布には 指数 分布 を選択する. 指数分布からのサンプリングには逆関数法が使える. : の一様乱数 s p (s) S = μ ​ e t −μ ​ s t ξ [0, 1] s = − ​ μ ​ t log(1 − ξ)

VREの評価( ) 距離 をサンプリングできたら, VREの各項を評価していく. Question 散乱光 をどう評価するか? レンダリング方程式で入射光 を評価するのと同様に, に対して再帰的にVREを適用して評価 発光と散乱のどちらを評価するか? ロシアンルーレットによってどちらを評価するか選択する(吸収/散乱イベントサンプリング) Term1 s Term1(s) = T(x, x − s )(μ ​ (x − ω a s )L ​ (x − ω e s , ) + ω ω μ ​ (x − s s )L ​ (x − ω s s , )) ω ω L ​ (x − s s , ) ω ω → L (x, ) i ω L ​ s →

イベントサンプリング ロシアンルーレットを用いることで, は以下のように評価できる. : 吸収イベントを選択する確率 : 散乱イベントを選択する確率 , は以下のように設定することで重点的サンプリングが行える. Term1(s) P ​ a P ​ s (s) = Term1 T(x, x − s )( ​ 1 ​ + ω P ​ a μ ​ L ​ (x − s , ) a e ω ω {ξ

散乱光の評価 : 位相関数. BRDFのボリューム版. モンテカルロ積分を適用すると 位相関数から方向 を重点的サンプリングして の方向にレイを飛ばし, の評価に再びVREを使う. L ​ (x, ) = s ω ​ f ​ ( , )L ​ (x, )dσ( ) ∫ S2 p ω ω′ i ω′ ω′ f ​ ( , ) p ω ω′ ​ (x, ) = L ​ s ω ​ p ​ ( ​ ) D ωi ′ f ​ ( , ​ )L ​ (x, ​ ) p ω ωi ′ i ωi ′ ω ω L (x, ​ ) i ω ​ i ′

位相関数 位相関数には Henyey-Greenstein phase function(HG 位相関数) がよく使われる. : と のなす角 : 散乱方向を制御するパラメーター HG位相関数からは解析的に方向 の重点的サンプリ ングが行える. θ ω ω′ g f ​ ( , ) = p ω ω′ ​ ​ 4π 1 (1 + g + 2g cosθ) 2 ​ 2 3 1 − g2 ω

VREの評価( ) は再びVREを適用するか, あるいは物体表面の場合はレンダリング方程式を評価すれば良い. まとめるとVREの評価は以下のようになる. : 距離サンプル : の一様乱数 Term2 Term2 = T(x, x − t )L(x − ω t , ) ω ω L(x − t , ) ω ω s ξ [0, 1] ​ ​ (x, ) L ω = ​ 1 ​ p ​ (s) S T(x, x − s )( ​ 1 ​ + ​ 1 ​ ) ω P ​ a μ ​ L ​ (x−s , ) a e ω ω {ξ

VREの評価 アルゴリズム的には以下のようにまとめられる. 1. 距離 をサンプリング 2. の場合 1. の一様乱数 を生成 2. なら を評価して終了 3. なら位相関数から を生成し、レイを まで進めて新たな方向 をセットし、1に戻る 3. の場合は を評価する。レイを まで進め、別の媒質に入る場合は1に戻る。真空 の場合は通常のパストレーシングで計算する s s < t [0, 1] ξ ξ < P ​ a L ​ (x − e s , ) ω ω ξ ≥ P ​ a ω′ x − sω ω′ s ≥ t L(x − t , ) ω ω x − tω

実装例 Mediumクラスで自由行程サンプリング, 散乱方向のサンプリングを行う. class Medium { protected: const std::shared_ptr phaseFunction; // 位相関数 public: Medium(float g) : phaseFunction(std::make_shared(g)) {} // 自由行程サンプリング, 散乱方向のサンプリングを行う. s > tの時はfalseを返す virtual bool sampleMedium(const Ray& ray, float distToSurface, Sampler& sampler, Vec3f& pos, Vec3f& dir, Vec3f& throughput) const = 0; // p1とp2の間の透過率を評価する virtual Vec3f transmittance(const Vec3f& p1, const Vec3f& p2, Sampler& sampler) const = 0; // 位相関数を評価する Vec3f evalPhaseFunction(const Vec3f& wo, const Vec3f& wi) const { return phaseFunction->evaluate(wo, wi); } };

実装例 波長非依存の一様媒質で発光が無い場合は sampleMedium は以下のように実装出来る. ` ` bool sampleMedium(const Ray& ray, float distToSurface, Sampler& sampler, Vec3f& pos, Vec3f& dir, Vec3f& throughput) const override { // 自由行程サンプリング const float t = -std::log(std::max(1.0f - sampler.getNext1D(), 0.0f)) / sigma_t; // s > tの場合 if (t > distToSurface - RAY_EPS) { pos = ray(distToSurface); dir = ray.direction; throughput = Vec3f(1); return false; } // 散乱方向をサンプリング phaseFunction->sampleDirection(-ray.direction, sampler, dir); pos = ray(t); throughput = Vec3f(sigma_s / sigma_t);

実装例 RayにMediumをスタックで保持させる. Rayが媒質に入るたびに mediums にMediumをpushする. 逆に媒質から出ていく時は mediums からpopする. class Ray { public: Vec3f origin; Vec3f direction; Vec3f throughput; std::stack mediums; }; ` ` ` `

実装例 純粋パストレーシングによるVREのモンテカルロ積分は以下のように実装できる. while (depth < maxDepth) { if (scene.intersect(ray, info)) { // ロシアンルーレット(省略) // 媒質がある場合は自由行程サンプリング, 散乱方向のサンプリングを行う bool is_scattered = false; if (ray.hasMedium()) { const Medium* medium = ray.getCurrentMedium(); Vec3f pos; Vec3f dir; Vec3f throughput_medium; is_scattered = medium->sampleMedium(ray, info.t, sampler, pos, dir, throughput_medium); // レイの更新 ray.origin = pos; ray.direction = dir; // throughputの更新

bool is_reflected_or_refracted = false; // s > tの場合 if (!is_scattered) { // 媒質終端が光源の場合 if (info.hitPrimitive->hasAreaLight()) { radiance += ray.throughput * info.hitPrimitive->Le(info.surfaceInfo, -ray.direction); break; } // 媒質終端が物体表面の場合は, BRDFから方向サンプリング Vec3f dir = ray.direction; if (info.hitPrimitive->hasSurface()) { float pdf_dir; const Vec3f f = info.hitPrimitive->sampleBxDF( -ray.direction, info.surfaceInfo, TransportDirection::FROM_CAMERA, sampler, dir, pdf_dir); // throughputの更新 ray.throughput *= f * cosTerm(-ray.direction, dir, info.surfaceInfo, TransportDirection::FROM_CAMERA) / pdf_dir; is_reflected_or_refracted = true; }

// レイが保持する媒質の更新 updateMedium(ray, dir, info); // レイの更新 ray.origin = info.surfaceInfo.position; ray.direction = dir; } // depthの更新 if (is_scattered || is_reflected_or_refracted) { depth++; } } else { // レイがシーン外に飛んでいった場合 break; } }

吸収/散乱係数が波長依存の場合 世の中に存在する多くの媒質は波長ごとに吸収/散乱係数が異なる. 水の吸収係数 Question 波長依存の場合, 自由行程サンプリング, 散乱方向サンプリングをどのように行えば良いのか? 波長を最初にサンプリングし, その後に自由行程サンプリング, 散乱方向サンプリングを行う →

波長サンプリング 個の波長の中から一様にサンプリングする方法はfireflyが出やすい. そこで, Rayのthroughputに比例するように波長をサンプリングする. 一様サンプリングとthroughput比例の比較 C

VREの評価 1つの波長だけでVREを評価すると, 他の波長が評価されないため色ノイズが顕著に出る. そこで, MISを用いて全波長分の寄与を一気に計算することを考える(Hero wavelength sampling). 1つの波長だけで評価する場合とMISの比較

Hero wavelength sampling 1つの波長 (Hero wavelength)を選んでパスを構築する. この時, 同様のパスを構築出来る戦略は波長の個数分 だけ存在する. 全波長分の寄与はMIS(バランスヒューリスティック)を用いて以下で計算できる. : X等色関数 : 波長 を選ぶ確率 : 波長 を使ってパス を構築する時の条件付き確率密度 λh (λ ​ ) x ˉ j P(λ ) h λh p( ∣λ ​ ) x ˉ k λ ​ k x ˉ ​ X = ​ ​ (λ ) j=1 ∑ C ​ P(λ )p( ∣λ ) ∑k=1 C k x ˉ k L( , λ ) x ˉ j x ˉ j

実装例 波長依存の一様媒質で発光がない場合は`sampleMedium`は以下のように実装できる. bool sampleMedium(const Ray& ray, float distToSurface, Sampler& sampler, Vec3f& pos, Vec3f& dir, Vec3f& throughput) const override { // 波長サンプリング Vec3f pmf_wavelength; const uint32_t channel = sampleWavelength( ray.throughput, (sigma_s / sigma_t), sampler, pmf_wavelength); // 自由行程サンプリング const float t = -std::log(std::max(1.0f - sampler.getNext1D(), 0.0f)) / sigma_t[channel]; // s > tの場合 if (t > distToSurface - RAY_EPS) { pos = ray(distToSurface); dir = ray.direction; const Vec3f tr = analyticTransmittance(distToSurface, sigma_t); const Vec3f p_surface = tr; const Vec3f pdf = pmf_wavelength * p_surface; throughput = tr / (pdf[0] + pdf[1] + pdf[2]);

// 散乱方向のサンプリング phaseFunction->sampleDirection(-ray.direction, sampler, dir); pos = ray(t); const Vec3f tr = analyticTransmittance(t, sigma_t); const Vec3f pdf_distance = sigma_t * tr; const Vec3f pdf = pmf_wavelength * pdf_distance; throughput = (tr * sigma_s) / (pdf[0] + pdf[1] + pdf[2]); return true; }

吸収/散乱係数の設定 吸収/散乱係数をアーティストが直接指定して見た目を制御するのは難しい. そこで, Disneyによって (Base Color)と (Scattering Distance)によって吸収/散乱係数を設定する方法が考 案されている. Blenderでも同様の方法が取られている. A d ​ ​ ​ μ ​ t μ ​ s s μ ​ t = 1 − exp −5.09406A + 2.61188A − 4.31805A 2 3 = 1.9 − A + 3.5(A − 0.8)2 = ​ ds 1

References Wilkie, Alexander, et al. "Hero wavelength spectral sampling." Computer Graphics Forum. Vol. 33. No. 4. 2014. Chiang, Matt Jen-Yuan, Peter Kutz, and Brent Burley. "Practical and controllable subsurface scattering for production path tracing." ACM SIGGRAPH 2016 Talks (2016): 1-2. Fong, Julian, et al. "Production volume rendering: Siggraph 2017 course." ACM SIGGRAPH 2017 Courses. 2017. 1-79. Wrenninge, Magnus, Ryusuke Villemin, and Christophe Hery. Path traced subsurface scattering using anisotropic phase functions and non-exponential free flights. Tech. Rep. 17-07, Pixar. https://graphics. pixar. com/library/PathTracedSubsurface, 2017. Novák, Jan, et al. "Monte Carlo methods for volumetric light transport simulation." Computer Graphics Forum. Vol. 37. No. 2. 2018.