解析的環境遮蔽

解析的に環境遮蔽を求めたい。


環境遮蔽(Ambient occlusion)を計算する場合、通常はモンテカルロ法を使った数値計算問題として解く。モンテカルロ法は解析解がとても複雑であったり、そもそも存在しないような場合であっても使えるためとても便利であるが、一方で標本の数が無限にならない限り真値に対する誤差は0にはならない。そもそも環境遮蔽の解析解が存在するならばモンテカルロ法を使う理由が弱くなる。いくつかの場合について解析解を求めてみる。

環境遮蔽

環境遮蔽は、表面上のある点$\bar{p}$がどの程度遮蔽されているかを表したものだ。具体的には半球上すべてが開けている状態をcosine項を考慮した上で1単位とする。cosine項を半球上に掛けると$\pi$になるので、遮蔽をcosine項で積分したのちに$\pi$で正規化すれば良いことになる。


念のため。

$ \begin{align} \int_{半球} cos \theta dA & = \int_{0}^{2\pi} \int_{\pi/2}^{0} cos\theta~sin\theta~d\theta d\phi \nonumber \\\ & = \int_{0}^{2\pi} \int_{\pi/2}^{0} \cfrac{1}{2} cos(2\theta)~d\theta d\phi \nonumber \\\ & = 2 \pi~\cfrac{1}{2} \nonumber \nonumber \\\ & = \pi \nonumber \\\ \end{align} $

また半球をcosineで射影した先は円になるので直感的にも$\pi$になることがわかる。


結局、環境遮蔽$A_{\bar{p}}$は次のようになる。

$ A_{\bar{\boldsymbol p}} = \cfrac{1}{\pi} \int_{半球} V(\hat{\boldsymbol \omega})~cos\theta~dA $

$V(\hat{\boldsymbol \omega})$は$\boldsymbol \omega$方向への可視性(visibility)で、遮蔽されていなければ1、遮蔽されていれば0を取る。この形は、計算したい点側から見た式で、モンテカルロ積分を行う場合は好都合だが、解析解を求めるときには使いづらい。解析解を求めるときは遮蔽する物体側から見た形の式の方が便利だ。まず微小面$dA$の${\bar p}$に対しての一般的な立体角は、そこへ向いたベクトル$\boldsymbol{p}$と、その微小面の法線$\boldsymbol{n_s}$から次のようになる。

$ \Omega = \int_{遮蔽物上} \cfrac{ \hat{\boldsymbol p} \cdot {\boldsymbol n_{s}} }{ ||{\boldsymbol p}||^2 } dA $

これに正規化項とcosine項を追加し、積分範囲を半球に制限すれば遮蔽する物体から見た形の環境遮蔽の式となる。

$ A_{\bar{p}} = \cfrac{1}{\pi} \int_{遮蔽物上} \cfrac{ \hat{\boldsymbol p} \cdot {\boldsymbol n_s} }{ ||{\boldsymbol p}||^2 } (\hat{\boldsymbol p } \cdot {\boldsymbol n})~dA $

これで環境遮蔽を解析的に求める下準備が出来た。

まず球体について考える。球体の中心を$P{s}$、半径を$R{s}$とする。$\bar{p}$から見て球が遮蔽している領域は円になる。つまり$A_{\bar{p}}$を求めるにはこの円上を積分すればよいことになる。

sphere1

この円の中心へのベクトルを$\boldsymbol{ P{c} }$、半径を$R{c}$とする。球に対する接線が球の中心から接線への半径方向と垂直で、$\boldsymbol P_{s}$と円の半径方向が垂直であることを使って次のように求まる。

$ R_{c} = \cfrac{R_{s}}{\|\boldsymbol P_{s}\|} \sqrt{\|{\boldsymbol P_{s}}\|^2 - R_{s}^2} \\\ $
$ \begin{align} \|{\boldsymbol P_{c}}\| & = \cfrac{R_c}{R_s} \sqrt{ {\|\boldsymbol P_{s}\|}^2 - R_{s}^2 } \nonumber \\\ & = \cfrac{1}{\|\boldsymbol P_{s}\|} (\|{\boldsymbol P_{s}}\|^2 - R_{s}^2) \nonumber \\\ \end{align} $

sphere2

円の極座標$(r,\phi)$を使って微小面積$dA$は次のように表現できる。

$ dA = r~d \phi~dr $

ここでもう一度$A_{\bar{p}}$がどういう形をしていたかを見てみる。

$ A_{\bar{p}} = \cfrac{1}{\pi} \int_{遮蔽物上} \cfrac{ \hat{\boldsymbol p} \cdot {\boldsymbol n_s} }{ ||{\boldsymbol p}||^2 } (\hat{\boldsymbol p } \cdot {\boldsymbol n})~dA $

このうち現時点で次の要素が求められる。

$ ||{\boldsymbol p}||^2 = { ||{\boldsymbol P_c}||^2 + r^2 } \\\ \hat{\boldsymbol p} \cdot {\boldsymbol n_s} = \cfrac{ ||{\boldsymbol P_c}|| }{ \sqrt{ ||{\boldsymbol P_c}||^2 + r^2 } } \\\ $

残りの$\hat{\boldsymbol p } \cdot {\boldsymbol n}$は次の図のようにして求める。

![sphere3]

結局$\hat{\boldsymbol p } \cdot {\boldsymbol n}$は次のようになる。

$ \begin{align} \hat{\boldsymbol p } \cdot {\boldsymbol n} & = \cfrac{ \boldsymbol p \cdot {\boldsymbol n} }{ ||\boldsymbol p|| } \\\ & = \cfrac{ ||P_c||(\hat{\boldsymbol P_c} \cdot {\boldsymbol n}) - r~cos(\phi)~sin(\alpha) }{ ||\boldsymbol p|| } \\\ \end{align} $

この結果を入れてみると次のようになる。

$ A_{\bar{p}} = \cfrac{1}{\pi} \int_{0}^{R_c} \int_{0}^{2\pi} \cfrac{ \hat{\boldsymbol p} \cdot {\boldsymbol n_s} }{ ||{\boldsymbol p}||^{5/2} } (||P_c||(\hat{\boldsymbol P_c} \cdot {\boldsymbol n}) - R_c~cos(\phi)~sin(\alpha))~r~d \phi~dr $

第二項目に注目する。$cos(\phi)$は$\pi$を原点にすると奇関数になり、他の部分は遇関数になっている。そのため第二項目は丸々キャンセルされる。残った第一項目を展開していく。

$ \begin{align} A_{\bar{p}} & = \cfrac{1}{\pi} \int_{0}^{R_c} \int_{0}^{2\pi} \cfrac{ \hat{\boldsymbol p} \cdot {\boldsymbol n_s} }{ ||{\boldsymbol p}||^{5/2} } (||P_c||(\hat{\boldsymbol P_c} \cdot {\boldsymbol n}) )~r~d \phi~dr \nonumber \\\ & = \cfrac{1}{\pi} \int_{0}^{R_c} \int_{0}^{2\pi} \cfrac{ ||\boldsymbol P_c|| }{ ({ ||{\boldsymbol P_c}||^2 + r^2 })^2 } (||P_c||(\hat{\boldsymbol P_c} \cdot {\boldsymbol n}) )~r~d \phi~dr \nonumber \\\ & = \cfrac{||\boldsymbol P_c||^2 (\hat{\boldsymbol P_c} \cdot {\boldsymbol n}) }{\pi} \int_{0}^{R_c} \int_{0}^{2\pi} \cfrac{ r }{ ({ ||{\boldsymbol P_c}||^2 + r^2 })^2 } ~d \phi~dr \nonumber \\\ & = \cfrac{||\boldsymbol P_c||^2 (\hat{\boldsymbol P_c} \cdot {\boldsymbol n}) }{\pi} 2\pi \int_{0}^{R_c} \cfrac{ r }{ ({ ||{\boldsymbol P_c}||^2 + r^2 })^2 } ~dr \nonumber \\\ & = 2||\boldsymbol P_c||^2 (\hat{\boldsymbol P_c} \cdot {\boldsymbol n}) \int_{0}^{R_c} \cfrac{ r }{ ({ ||{\boldsymbol P_c}||^2 + r^2 })^2 } ~dr \nonumber \\\ & = -||\boldsymbol P_c||^2 (\hat{\boldsymbol P_c} \cdot {\boldsymbol n}) \Biggl[ \cfrac{ 1 }{ { ||{\boldsymbol P_c}||^2 + r^2 } } \Biggr] ^{R_c}_0 \nonumber \\\ & = -||\boldsymbol P_c||^2 (\hat{\boldsymbol P_c} \cdot {\boldsymbol n}) \Biggl( \cfrac{ 1 }{ { ||{\boldsymbol P_c}||^2 + R_c^2 } } - \cfrac{ 1 }{ { ||{\boldsymbol P_c}||^2 } }\Biggr) \nonumber \\\ & = (\hat{\boldsymbol P_c} \cdot {\boldsymbol n}) \cfrac{ ||{\boldsymbol P_s}||^2 - R_s^2 }{ ||{\boldsymbol P_s}||^2 } \cfrac{ R_s^2 }{ { ||{\boldsymbol P_s}||^2 - R_s^2 } } \nonumber \\\ & = (\hat{\boldsymbol P_s} \cdot {\boldsymbol n}) \Biggl( \cfrac{R_s}{||{\boldsymbol P_s}||} \Biggr)^2 \nonumber \\\ \end{align} $

これで球の環境遮蔽が解析的に求まった。

$ A_{\bar{p}} = (\hat{\boldsymbol P_s} \cdot {\boldsymbol n}) \Biggl( \cfrac{R_s}{||{\boldsymbol P_s}||} \Biggr)^2 $

注意点としてはこの計算は${\hat p}$の半球上から球全体が見えていることを仮定していることだ。一部でも半球上から隠れてしまった場合正しい結果にはならない。$(\hat{\boldsymbol Ps} \cdot {\boldsymbol n})$の部分からもわかるように、ちょうど半分隠れたときに$A{\hat p}$は0になってしまうが、現実には球全体が隠れきった時点で0になるのが正しい。同様に複数球を並べ、それらの環境遮蔽を単純に足しあげる方法も他の球に一部でも遮蔽されている方向があれば正確ではなくなる。全く遮蔽されない状況は平面に対して平行な面上に球が重なることなく分布しているときなど、かなり特殊な状況なので実際上は複数個置いた時点で正確な結果は期待できない

n-gon

次はn-gonについて考える。$A_{\bar{p}}$の式を再掲する。

$ A_{\bar{p}} = \cfrac{1}{\pi} \int_{遮蔽物上} \cfrac{ \hat{\boldsymbol p} \cdot {\boldsymbol n_s} }{ ||{\boldsymbol p}||^2 } (\hat{\boldsymbol p } \cdot {\boldsymbol n})~dA $

この式は「半球上で面積分している」わけだが、Stokesの定理から考えれば、この積分の中の式が回転になる関数をうまいこと見つけだせれば経路積分に変形することが出来る。


Stokesの定理。「ある式の回転の面積分」は「ある式の閉経路積分」となる定理。

$ \int \int_{面} \nabla \times {\boldsymbol F} \cdot d{\boldsymbol S} = \oint_{囲む経路} {\boldsymbol F} \cdot d{\boldsymbol l} $

n-gonに関する場合のその式がBaumによって発見された[1](導出は[2])。それを使い再度$A_{\hat p}$を書きなおすと次のようになる(本来はラジオシティのForm-Factorとして発見されたものだが話の繋がり上、環境遮蔽の式として書いた)。

$ A_{\hat p} = \cfrac{1}{2 \pi} \Sigma_{全エッジ} {\boldsymbol N} \cdot \cfrac{ {\boldsymbol R_i} \times {\boldsymbol R_{i+1}}}{||{\boldsymbol R_i} \times {\boldsymbol R_{i+1}}||} \gamma_i $

$\boldsymbol R_i$はi番目の頂点座標ベクトル、$\gamma_i$はi番目とi+1番目頂点がなす角度だ。$N$と$\hat {\boldsymbol Vi}$と$\hat{\boldsymbol V{i+1} }$が半球上に張る領域に相当し、それらを加減算(頂点の順番の回転方向が逆になった場合に減算になる)したものと考えると、そこそこ感覚的に理解できる(気がする)。またn-gonを構成する三角形が同一平面内にあることを仮定していないので、同一平面内になくてもこの式は成立する。

凸包

n-gonに関する環境遮蔽が計算できたので、任意のメッシュに対する環境遮蔽も計算できるような気がしてくるが、それは正しくない。なぜならあるn-gonが他のn-gonを遮蔽するようなシチュエーション、つまりダブルカウントをしてしまう場合を考えないといけないからだ。凸包であればあるn-gonが他のn-gonを遮蔽している場合、遮蔽される側のn-gonが必ず裏を向いているため、この制限を回避できる。凸包の環境遮蔽を計算するには次のような手順を踏む。

  1. 凸包中のあるポリゴンが${\hat p}$に対して表裏のどちらを向いているかを判定する。これは$\boldsymbol V_0$,$\boldsymbol V_1$,$\boldsymbol V_2$の作る面の法線と、${\boldsymbol V_0} - \hat{p}$の内積の正負から判定が出来る。条件分岐をさせずに結果を階段関数に掛けて後段は投機的に実行した方が早い場合もあるかもしれない。
  2. もし表を向いていればそのポリゴンの環境遮蔽を求める
  3. 2の結果を足し合わせる。

平面

最後に平面の環境遮蔽を算出する。この形状は積分を行わずにとても簡単に感覚的に算出することができる。この平面の法線を$\boldsymbol N_p$とする。平面は無限に大きい四角形なので、半球上で遮蔽されていると判定される領域の境界線とその平面は平行になる。

plane0

これを俯瞰すると次のような図になる。

plane1

この図の右半分にある楕円の面積は長径と短径の積に比例し、短径は${\boldsymbol N_p} \cdot {\boldsymbol N}$になることがわかる。まとめると平面による環境遮蔽は次のようになる。

$ A_{\hat p} = \cfrac{ {\boldsymbol N_p} \cdot {\boldsymbol N}+1}{2} $

結果

中央のスライダーをドラッグして比較することが出来ます。

球+16spp


## 正十二面体+16spp


参考

  1. http://pdf.aminer.org/000/593/020/improving_radiosity_solutions_through_the_use_of_analytically_determined_form.pdf
  2. http://www.multires.caltech.edu/pubs/fftr.pdf
  3. http://ambientocclusion.hatenablog.com/entry/2013/10/15/223302
  4. http://iquilezles.org/www/articles/sphereao/sphereao.htm
  5. http://en.wikipedia.org/wiki/Solid_angle
  6. https://www.shadertoy.com/view/XdjSDy
  7. https://www.shadertoy.com/view/4djXDy
  8. https://www.shadertoy.com/view/4djSDy

[sphere3]: /img/analytic_ambient_occlusion/sphere3.png