VC/GCAD 合同シンポジウム 2015

Narrow beam 理論を用いた散乱媒質中の物体の画像生成
Rendering Objects in Random Scatering Media with the Narrow Beam Theory
1)4)

2)4)

1)

1) 東邦大学
2) 北海道大学
3) 修道大学
4) UEI リサーチ

Toho University
Hokkaido University
Shudo University
UEI Research

3)4)

新谷 幹夫
土橋 宜典
白石 路雄
西田 友是
Mikio SHINYA, Yoshinori DOBASHI, Michio SHIRAISHI, Tomoyuki NISHITA

1

はじめに散乱媒質中で光は粒子により散乱し拡がっていく。これは、視線が散乱により拡がるとも解釈することができる (図金属や結晶などを除くほとんどの物体は散乱性の反射特は定式化も容易で、古くから depth-cue の一種としても用

2)。物体反射光の散乱を扱うためには、視線の散乱による拡がりを評価し、これに基づき「にじみ」効果を表現するのが合理的である。このためには、散乱媒質内での光束の拡がりを陽に求める必要がある。光束の拡がりは基本的な量であるにもかかわらず、これまで CG 分野では光束の拡がりを効率的かつ明示的に求める手法はなかった。
そこで本研究では、Narrow beam 理論 [2] を画像生成に導入し、物体反射光の散乱を直接的かつ効率的に計算する手法を提案する。Narrow beam 理論は、光強度は１方向にいられている。照明光の散乱は、いわゆる表面下散乱であ集中していることを前提とした散乱近似理論である。視線り、様々な手法が開発されている。ところが、物体光の散の散乱を扱う場合、もとの視線方向に集中しているので、乱にフォーカスした研究例をあまりみない。path-tracing

適しているといえる。理論を基に、光線入射に対する拡散性を持ち、多重散乱の表現は写実的画像生成の重要な課題
1 物である。散乱現象における見かけ上の特徴としては、°2 照明光の散乱（にごり）
3 物体体光の減衰（かすみ）、°、°光の散乱（にじみ）、などがあげられる。例えば、透明度の2 ）低い池を考えると、水は白く濁って見え（°、池の石や魚1 ）は深さに応じてコントラストが低下し（°、境界や模様3 ）がボケて見える（°。このうち、物体からの反射光の減衰[1] などの確率的手法を用いれば表現が可能であるが、効率的かつ有効な手法が必要とされている。
散乱光の画像生成は、光輸送方程式を解くことに帰着され、その解法により手法を分類することができる。統計的アプローチでは、光線を確率的にサンプルし、モンテカルロ積分により解を求める。path-tracing などが代表例であ式を導出した。ついで、画像生成への応用に適するように、散乱分布をパラメータ表現する定式化を考案した。これにより、closed-form の単純な計算で、奥行値、散乱パラメータから拡がり関数を算出できる。フォトントレーシング結果と提案手法との比較実験を行ったところ、よい一致を見た。この定式化を用いて、物体散乱光の計算を画像空間上る。充分なサンプルを取れば収束するので、物体光の散乱で行うアルゴリズムを提案し、GPU 実装し、予備的な画も産出されるが、一般に処理時間を要する。拡散近似アプ像生成実験を行った。本研究の貢献は以下のようである。

ローチは、散乱光の方向性が弱いことを前提に導かれる拡散方程式を解くことで散乱場を求める。dipole/multi-pole 法• 物体光の散乱を効率的に計算する手法を提案した[3, 4] やその拡張手法 [5, 6] は、対象媒質を層状であると仮定し、拡散方程式の解析解を基に BSSRDF(Bi-directional• Narrow beam 理論により光束の散乱媒質中の拡がりが求められることを示した

Sub-Surface Reflectance Distribution Function) を求め、照明光の表面下散乱を表現する。物体反射光を２次光源と• Narrow beam 理論を基に、画像生成で利用しやすい解析的近似関数を提案した

して扱えば、その散乱も扱えるが、視線による物体像の変化が表現できず、見通しは悪い。また、層状仮定に反する場合には合理性が低く、精度も期待できない。ray-marching

2アプローチでは、視線上の経路積分で視線方向の散乱光を求める。１次散乱を求める場合の定石であるが、多重散乱物体光の散乱
図 1 に示すように、散乱媒質中の物体の像はぼやけて見も局所的多重反射 [7] や plane-parallel 近似 [8] などと組み1 物体表面での反射を光輸送える。この計算のためには、°合わせることで計算可能である。しかしこのアプローチに2 物体反射光を２次過程に組み込み、統計的手法を適用、°おいても、視線は直線であるので、物体光の散乱を扱うこ光源として扱い、multipole 法 [4] など拡散近似手法を適とができない。

3 視線の散乱による物体表面近傍での拡がりを求め、用、°13:1

VC/GCAD 合同シンポジウム 2015

0scattering m edium‫ݖ‬Ƹz‫ݎ‬Ƹrߴs‫ݏ‬Ƹ

図 3: 座標系と実験配置
図 1: 物体光の散乱. 栗がボケて見える.
表 1: Symbols.

BC㼐㼕㼜㼛㼘㼑㻌r̂,r, z

位置 r̂ = r + z ẑŝ,s方向 ŝ = s + ẑσt

散乱係数σs

散乱係数σa吸収係数

ρn

粒子密度σs = σt − σs-κA㼐㼕㼜㼛㼘㼑㻌

r に対するqs に対するフーリエ変数+フーリエ変数

p(s)

位相関数

P (κ)

p のフーリエ変換

I0 (r, s)

入射光分布F0 (κ, q)

I0 のフーリエ変換

図 2: 視線の散乱と multi-pole モデル

Z
物体反射光をサンプル、などのアプローチが考えられる。

p(ŝ, ŝ0 )I(r̂, ŝ0 )dŝ0 (1)ρn (σs /4π)

1 は自然に物体光散乱を表現するが、処理量がアプローチ°4π2大きいのでより効率的な手法が期待される。アプローチ°と表す。ここで、ŝ は方向、r̂ は位置を表す３次元ベクトは、拡散近似の性格上、方向性の表現能力が乏しく、空間ル、ρn は密度、σt , σs は減衰係数および散乱係数、p は位構造が重要な物体像の生成には不向きである。図 2 に示

相関数で、積分範囲は単位円上である。I を直接光成分 Iri

すように、散乱媒質中の物体点 A を２次光源と考えると、と間接光成分 Id に分解し、multi-pole 法では物体法線方向に仮想光源が置かれること

I = Iri + Idになり、媒体表面の点 B が A 点の像の中心となる。本来、とする。以下の２点A 点方向の視線上の点 C に A 点の像の中心が期待されるが、物体の像が媒質表面に投影されて貼りつき、媒質と物• 光強度は z 軸方向に集中していると仮定する。負方向の伝搬は無視する。（ŝ · ẑ = cos θ ' 1）

体との視差が反映されないことになる。肌のように視差が無視できる場合以外は大きな問題である。この点は、改良• 位相関数の変化が散乱強度の高い範囲では小さく、２次近似できる。（式 1 の位相関数との積分で Id を２

手法 [5, 6] でも同様である。
3 の立場をとる。このために本研究では、アプローチ°次近似できる）

は、散乱による視線の拡がりを効率的に評価する必要があるが、Narrow beam 理論を適用することで実現できるこを仮定し、とを示す。

r̂ = r + z ẑ, ŝ ' s + ẑとする。ただし、ẑ は z 軸方向の単位ベクトルである。こ

3れにより、方向に関する球面積分を平面積分に置き換える。

Narrow beam 理論

また、散乱は等方的であるとし、本研究で導入する Narrow beam 理論の概略を示す。なp(s, s0 ) = p(|s − s0 |)お、使用する主なシンボルを表 1 に示す。支配方程式であとする。これにより、積分偏微分方程式が２階線形偏微分る光輸送方程式を(∇ · ŝ)I(r̂, ŝ) =(2)方程式に簡略化され、フーリエ変換などを用いて、計算を−ρn σt I(r̂, ŝ) +進めると、13:2

VC/GCAD 合同シンポジウム 2015

Iri (z, r, s)=(1/16π 4 )

Z ∞

P (q)

Z ∞
dκ

−∞

exp(−iκ · r)

=4π exp(−Bq 2 )(10)とする。ここで、z 軸方向の光線の入射は−∞exp(−is · q)F0 (κ, q + κz) exp(−σt ζ(z))dqZ
I0 (r, s) = δ(r)δ(s)Id (z, r, s) = (1/2π)4 exp(−ρn σa )z Fd (z, κ, q + κz)とデルタ関数により表せる。すると、Zexp(−is · q)dq exp(−iκ · r)dκ(3)F0 (κ, q) = 1(11)Z z00Fd (z, κ, q) = (σt /4π)F0 (κ, q)
ρn (z )P (q − z κ)0となる。これを式 3 に代入し、整理すると線積分exp(−σs ζ(z 0 )) exp(−q 2 Â(z, z 0 ))dz 0(4)Z z
Id (z, r, s) = KG (z) exp{−s2 /(4Σ2G (z))}
Z zζ(z) =
ρn (z 0 )σt dz 0(5)0(π/Σ2G ) exp(−ρn σs z 0 )dz 0 (π/C 2 (z 0 ))Z zÂ(z, z 0 )A(z 00 )dz 00=0(6)exp{−(r − zs + Bz 0 s/Σ2G )2 /(4C 2 (z 0 ))}(7)= KG (z) exp{−s2 /(4Σ2G )} · π 2Z zdz 0 exp(−ρn σs z 0 )/(Σ2G (z 0 ) · C 2 (z 0 )) ·z0= ρn (z)(σs /4) < θ2 >ZZ
22<θ > =θ p(s)ds/p(s)dsA(z)4π(8)

4π0exp{−(r + (−z + Bz 0 /Σ2G (z 0 )s)2 /(4C 2 (z 0 ))}と位置 (r, z)、方向 s の散乱光強度が求まる。ただし、i は虚数単位、F0 は入射光 I0 (r, s) のフーリエ変換
ZZ
F0 (κ, q) = exp(iq · s)ds I0 (r, s) exp(iκ · r)dr(12)を得る。ここで、P は位相関数 p のフーリエ変換である。また、κ,q は r,sに対応するフーリエ変数であり、これらに関する積分は平Σ2G (z 0 ) =20Bza =ρn σsC (z ) =面全体での面積分である。導出の詳細は文献 [2] を参照されたい1 。

KG (z) =4

視線拡がり関数

B + A(z − z 0 )02(13)[1 − B/Σ2G ]

(14)
(15)

4(1/2π) σs exp(−ρn σa z)(16)

また、直接光成分は、

Iri (z, r, s) = δ(r − zs)δ(s) exp(−ρn σt z)

式 3 はフーリエ変換と視線方向の線積分により直接計算できるが、画像生成応用を考えるとより簡便な定式化が望(17)となる。導出過程を補足資料 1 に示す。

まれる。本節では、まず、入力を光線（デルタ関数）、位相関数をガウス関数とすることで、フーリエ変換を解析的に行い、視線上の積分により直接計算することが可能であ4.2モーメントることを示す。

Id の 0 次モーメント φ、２次モーメント RZ Z
φ(z) =Id (z, r, s)dsdr
Z Z
R(z) =
ρ2 Id (z, r, s)dsdr

各画素毎に拡がり領域内のサンプル点に関して視線積分を行うのはできれば避けたい。そこで、散乱光の２次モーメントと０次モーメント（総量）を計算し、これらのパラメータから先験的関数を規定することを考える。濃度 ρn(18)(19)が一定値である場合には、モーメントは奥行値から解析的に計算できることを示す。先験的関数としては、１次散乱を求める。なお、Id は対称関数であるので１次モーメントによる散乱分布を参考に正規積分関数 (erf) を用いた関数は 0 である。式 12 を式 19 に入れ、整理すると、形を提案する。

4.1

光線入力

=exp(−ρn σa z)(1 − exp(−ρn σs z))

(20)

R(z)

=

4 exp(−ρn σa z) ·(21)22[exp(−ρn σs z){Az /a − 2B/a } +位相関数 p(s) がガウス関数、

p(s)

φ(z)

=

(1/B) exp(−s2 /4B)

{2B/a2 − 2Bz/a + (−A/a + B)z 2 + Az 3 }](9)を得る。導出過程を補足資料 2 に示す。

1 文献の式 (13-39) にタイポがあると思われるので注意を要する。

13:3

VC/GCAD 合同シンポジウム 2015

4.3

拡がり関数

計算されたモーメントを基に、先験的拡がり関数を決める。尚、導出過程は補足資料 3,4 に示す。
前節と同一条件で１次散乱光 I1 のみを計算すると、
Z z
I1 (z, r, s) = (σt /4π)exp(−σt ρn y)dy0
ZI0 (z − y, r + ys, s0 )p(s, s0 )ds0
(22)
4π
Z z
= (σt /4π) exp(−σt ρn z)
δ(r + ys)p(s)dy݀ߠ ൌ ͷι݀ߠ ൌ ͳͲι0となる。位置 (r, z) における全散乱光Z
Ψ1 (z, r) = I1 (z, r, s)ds(23)|r|݀ߠ ൌ ͵Ͳι

2.5を計算すると、正規積分関数√ Z ∞exp(u2 /2)duG(x) = 1/ 2π(24)を用い、

Ψ1 (z, r)1.05.0zx

estim atedphoton tracing=(σs ) exp(−σt ρn z) · (1/B)√√(1/|r|)(1/ π)G(|r|/(z 2B))݀ߠ ൌ ͸Ͳι

10-3(25)

図 4: 比較実験結果と表される。これを参考に、先験的な拡がり関数として、

RM S
0.7

Ψ(r) = α(1/|r|)G(|r|/β)(26)

0.60.5を採用する。

direct0.4

single scattering

パラメータ α,β はモーメント φ,R から0.3

proposedGauss

α = 2φ/β 2(27)β 2 = 3(R/φ)

(28)

0.2
0.1
001020と求めることができる。

3040506070

݀ߠ (degree)

図 5: 誤差と位相関数の拡がり角 dθ5比較実験くなることが示されている。photon-tracing 法の結果（左側）と提案手法の結果（右側）を比較すると、概ねよい一先験的拡がり関数の近似精度を確認するため、photon-致をみている。

tracing 法との比較実験を行った。図 3 に示すように、散図 5 は位相関数の拡がり角に対する root-mean-square

乱媒質は z の正領域にあり、原点から z 軸方向に光線を入(rms) 誤差を示している。誤差は 0.2 以下、概ね 0.1 前後であり、比較的よい近似が得られている。多重散乱の効果を射した。フォトンを確率的に散乱させ、通過した光線をカウントすることにより分布を求めた。一方、式 21,26 から見るため、直接光、１次散乱光成分を式 17 および 25 を用拡がり関数も求めて比較した。

いて計算し、rms 誤差を求めた。この結果も図 5 に示す。

σt = 1,σs = 0.95 に固定し、位相関数はガウス関数として、その拡がり角 dθ2 = 4B を変化させた。得られた分布

図のように提案手法では、多重散乱を考慮することにより、と計算された拡がり関数の比較を図 4 に疑似カラー表示で誤差が大きく減少していることが分かる。
なお、拡がり関数の関数系としてガウス関数示す。ただし、拡がり関数 Ψ(r)（式 26）は r = 0 で特異性があるので、その周回積分値 2πrΨ(r) を示している。深さΨrg (r) = α0 exp(−|r|2 /β 0 )z に従い、光線が拡がっていく様子が見られる。また、位相関数の拡がり角を大きくとると、散乱光の拡がりも大きなども試したが、提案関数が最良の結果を示した。参照の13:4

VC/GCAD 合同シンポジウム 2015

点を A とし、A を含む s1 に垂直な面を M とし、M と s2ʹߨ‫ݎ‬Ȳ0.5との交点を B とする。それぞれの視線が近傍画素に対応す

albedo= 0.70.45るとすれば、A,B における視線拡がり関数はほぼ等しいと0.40.35z=10.3(photon)

0.25z=20.2(photon)

仮定できる。拡がり関数を Ψ(r)、A における反射輝度をI 0 (A) とすれば、s2 の視線方向には、A の反射光が散乱により Ψ(|AB|)I 0 (A) だけ加算されることになる。したがって、s2 で観測される物体光の散乱強度 I(B) は、物体と交差する視線の物体点 A に関する和z=30.15(photon)

0.1

0.05
000.2

0.4

0.60.81|r|

I(B) =X

Ψ(|AB|)I 0 (A)(29)

Aʹߨ‫ݎ‬Ȳ

0.5

albedo= 0.99により与えられる。フィルタカーネル Ψ は、媒質表面から0.450.40.35z=10.3(photon)

0.25z=20.2(photon)の距離、散乱パラメータより式 20, 21, 26, 27, 28 を用いて算出される。
そこで、媒質表面の奥行画像、物体の奥行画像、物体の反射画像を入力として与えれば、式 29 による線形フィル

z=30.15(photon)0.1

タリングにより物体光の散乱光強度が算出される。これと0.05
0

既存手法による照明光の散乱強度、すなわち表面化散乱強|r|
00.2

0.4

0.6

0.81度とを合成し、画像を生成できる。前述のように、拡がり関数の計算は、式 21,26 を評価することのみで求まるので、

図 6: 誤差とアルベド処理コストは式 29 のフィルタリングのみであり、対話的速度で実行可能である。

ため、図 5 にガウス関数に対する rms 誤差を示す。

図 8 に生成画像例を示す。照明光の散乱は Plane-parallel

アルベド (σs /σt ) の影響を見るため、低アルベド媒質近似 [8] を用いて計算している。左側に入力として用いた

σs = 0.7, σt = 1 と高アルベド媒質 σs = 0.99 に対し、同様の比較を行った。位相関数の拡がり角を 30 度に固定し、z = 1, 2, 3 における計算結果を図 6 にプロットする。図に示すように、どちらのアルベドに対しても、photon-tracing法の結果とよく一致していることが分かる。

6

物体反射光と媒質の散乱光を示し、右側に提案手法による生成画像を示す。物体光のにじみが表現されていることが分かる。処理時間は、Windows PC (Intel Core i7-3770K@3.50GHz, GeForce GTX 690) 上で 170ms であった。

画像生成㼟㻞㻌㼟㻝㻌

Ȳ

図 8: 画像生成例

AB
M

7まとめ
本研究では、散乱媒質中の物体光の散乱に着目し、これ

図 7: 視線の拡がりとフィルタリングを画像空間で効率的に計算する手法を提案した。このため、視線の拡がり関数を用いた画像空間のアルゴリズムを構視線の拡がりを評価する理論として、Narrow beam 理論を築し、GPU 実装した。図 7 に示すように、視線 s1 は物体導入し、拡がりの算出方法を定式化した。また、画像生成に交差し、s2 が交差しない場合を考える。s1 と物体の交に適用するため、先験的な拡がり関数を導入し、視線の散13:5

VC/GCAD 合同シンポジウム 2015

乱分布を closed-form で記述することに成功した。photontracing 法との比較実験により、提案した拡がり関数が良好な精度を持つことを示した。
実装の完成度を高めて複雑な情景に適用するとともに、以下の拡張を行っていく。

• 媒質表面での屈折: 入力画像の生成に光線追跡法を適用することで表現可能• 不均質媒体への適応:視線毎の線積分で不均質媒体へも適応可能であり、ray-marching により実装可能• 微細凹凸のある媒質境界面への適用：表面からの入力光分布を加味することで実現可能
また、フィルタ処理が処理時間の大半を占めるので、多重解像度画像の利用や separable filter による近似など、効率的な実装も検討したい。
本稿では物体光の散乱を対象課題としたが、光束の散乱による拡がりは様々な局面で利用可能である。照明光に適用すれば、ヘッドライトや木漏れ日などの light shaft の多重散乱効果などの効率的表現が期待できる。また、散乱パラメータの測定への応用も検討したい。

参考文献[1] J. T. Kajiya, The rendering equation,
GRAPH’86, No. 4, pp. 143-150, 1986.

SIG-[2] A. Ishimaru, Wave propergation and scattering inrandom media, volume 1, Academic Press, NewYork, 1978.
[3] Jensen, H. W. and Marschner, S. R. and Levoy, M.
and Hanrahan, P., A practical model for subsurface
light transport, SIGGRAPH 2001,pp.511-518, 2001.
[4] Donner, C. and Jensen, H. W., Light diffusion in
multi-layered translucent materials, ACM Transactions on Graphics, vol. 24, No. 3, pp.1032-1039, 2005.
[5] E. d’Eon, D. Luebke, E. Enderton, Efficient Rendering of Human Skin, Proceedings of the Eurographics Symposium on Rendering Techniques, Grenoble,France, 2007
[6] E. d’Eon, G. Irving, A quantized-diffusion model for
rendering translucent materials, ACM Transactions
on Graphics, vol. 30, No. 4, 56, 2011.
[7] A. Zinke, C. Yuksel, A. Weber, and J. Keyser, Dual
scattering approximation for fast multiple scatteringin hair, ACM Transactions on Graphics vol. 26, No.
3, pp. 1-10 (2008).

13:6[8] M. Shinya, M. Shiraishi, Y. Dobashi, K. Iwasaki, T.
Nishita, A simplified plane-parallel scattering modeland its application to hair rendering, Pacific Graphics 2010, pp.85-92, 2010.