バリアンスシャドウマップのモーメント計算の行間を補う

はじめに

GPU Gem3 - 8.4 Variance Shadow Mapsで、シャドウマップのモーメントを計算している箇所があるんですが、
そこを行間を補って丁寧に解説するという記事です。

実際のアルゴリズム

頂点シェーダで深度値(光源から見た時のZ値)を計算しoutします。
フラグメントシェーダで補間された深度値がinされます。その値を $\mu$ とします。
現在のフラグメント周りでの深度値が $x,y$ によって決まるとします。
要するに $z=f(x,y)$ とします。
$xy$ 座標系の原点を現在のフラグメントとすれば、関数 $f$ は以下のように書けます。

$\begin{eqnarray} f(x,y)=\mu+\frac{\partial f}{\partial x}x+\frac{\partial f}{\partial y}y\tag{1} \end{eqnarray}$

式 $(1)$ の $\dfrac{\partial f}{\partial x},\dfrac{\partial f}{\partial y}$ はフラグメントシェーダが計算してくれるので心配ありません。
通常の画像処理なら、隣り合う画素から偏微分の値を計算できますが、今回は隣り合う画素の値を取得できない状態です。
ですが、同じことをもう一度言いますが、 $\dfrac{\partial f}{\partial x},\dfrac{\partial f}{\partial y}$ はフラグメントシェーダが計算してくれるので心配ありません。

チェビシェフの不等式を使うには、深度値の分散と平均が必要です。
その為に $x,y$ の確率分布を定めましょう。
計算しやすいように $x,y$ は相関が無い正規分布とします。
平均は原点 $(0,0)$ で標準偏差は $\sigma=1/2$ とします。分散は $\sigma^2=1/4$ です。
この時、 $x,y$ の確率分布の確率密度関数は以下のようになります。

$\begin{eqnarray} p(x,y)={\mathcal N}\left(\begin{pmatrix}x\\ y\end{pmatrix}\ \middle|\ \begin{pmatrix}0\\ 0\end{pmatrix},\ \begin{pmatrix}1/4 & 0\\ 0 & 1/4\end{pmatrix}\right)\tag{2} \end{eqnarray}$

式 $(2)$ は2次元正規分布ですが、 $x,y$ は無相関なので以下のように書けます。

$\begin{eqnarray} p(x,y)&=&p(x)p(y)\\ &=&{\mathcal N}\left(x\ \middle|\ 0,1/4\right){\mathcal N}\left(y\ \middle|\ 0,1/4\right)\tag{3} \end{eqnarray}$

これから先で期待値の計算をしますが、その為に簡単な計算を済ませておきます。
期待値の計算はうっとおしいかも知れませんが、 ${\rm E}[\cdot]_{p(x)}$ のように閉じ括弧の下側に確率分布を明示することにします。

$\begin{eqnarray} &&{\rm E}[\mu]_{p(x,y)}=\mu{\rm E}[1]_{p(x,y)}=\mu\tag{3}\\ &&{\rm E}[x]_{p(x,y)}={\rm E}[x]_{p(x)p(y)}={\rm E}[x]_{p(x)}=0\tag{4}\\ &&{\rm E}[y]_{p(x,y)}={\rm E}[y]_{p(x)p(y)}={\rm E}[y]_{p(y)}=0\tag{5}\\ &&{\rm E}[xy]_{p(x,y)}={\rm E}[xy]_{p(x)p(y)}={\rm E}[x]_{p(x)}{\rm E}[y]_{p(y)}=0\cdot 0=0\tag{6}\\ &&{\rm E}[x^2]_{p(x,y)}={\rm E}[x^2]_{p(x)p(y)}={\rm E}[x^2]_{p(x)}=\sigma^2+{\rm E}[x]_{p(x)}^2=1/4+0^2=1/4\tag{7}\\ &&{\rm E}[y^2]_{p(x,y)}={\rm E}[y^2]_{p(x)p(y)}={\rm E}[y^2]_{p(y)}=\sigma^2+{\rm E}[y]_{p(y)}^2=1/4+0^2=1/4\tag{8}\\ \end{eqnarray}$

(下の計算も丁寧にやっているので式 $(3)$ ～ $(8)$ は不要ですね。計算苦手な人のために残しておきます。)
さて、準備は終わりました。後は ${\rm M}_1,{\rm M}_2$ を計算すれば終わりです。

1次モーメント ${\rm M}_1={\rm E}[f]_{p(x,y)}$ を計算します。

$\begin{eqnarray} {\rm E}[f]_{p(x,y)}&=&{\rm E}\left[\mu+\frac{\partial f}{\partial x}x+\frac{\partial f}{\partial y}y\right]_{p(x,y)}\\ &=&{\rm E}\left[\mu+\frac{\partial f}{\partial x}x+\frac{\partial f}{\partial y}y\right]_{p(x)p(y)}\\ &=&{\rm E}\left[\mu\right]_{p(x)p(y)}+{\rm E}\left[\frac{\partial f}{\partial x}x\right]_{p(x)p(y)}+{\rm E}\left[\frac{\partial f}{\partial y}y\right]_{p(x)p(y)}\\ &=&\mu{\rm E}\left[1\right]_{p(x)p(y)}+\frac{\partial f}{\partial x}{\rm E}\left[x\right]_{p(x)}+\frac{\partial f}{\partial y}{\rm E}\left[y\right]_{p(y)}\\ &=&\mu\cdot 1 + \frac{\partial f}{\partial x}\cdot 0 +\frac{\partial f}{\partial y}\cdot 0\\ &=&\mu\tag{9} \end{eqnarray}$

2次モーメント ${\rm M}_2={\rm E}[f^2]_{p(x,y)}$ を計算します。

$\begin{eqnarray} {\rm E}[f^2]_{p(x,y)}&=&{\rm E}\left[\left(\mu+\frac{\partial f}{\partial x}x+\frac{\partial f}{\partial y}y\right)^2\right]_{p(x,y)}\\ &=&{\rm E}\left[\mu^2+\left(\frac{\partial f}{\partial x}\right)^2x^2+\left(\frac{\partial f}{\partial y}\right)^2y^2+2\mu\frac{\partial f}{\partial x}x+2\mu\frac{\partial f}{\partial y}y+2\frac{\partial f}{\partial x}x\frac{\partial f}{\partial y}y\right]_{p(x)p(y)}\\ &=&{\rm E}\left[\mu^2\right]_{p(x)p(y)}+{\rm E}\left[\left(\frac{\partial f}{\partial x}\right)^2x^2\right]_{p(x)p(y)}+{\rm E}\left[\left(\frac{\partial f}{\partial y}\right)^2y^2\right]_{p(x)p(y)}+{\rm E}\left[2\mu\frac{\partial f}{\partial y}y\right]_{p(x)p(y)}+{\rm E}\left[2\mu\frac{\partial f}{\partial y}y\right]_{p(x)p(y)}+{\rm E}\left[2\frac{\partial f}{\partial x}x\frac{\partial f}{\partial y}y\right]_{p(x)p(y)}\\ &=&\mu^2{\rm E}\left[1\right]_{p(x)p(y)}+\left(\frac{\partial f}{\partial x}\right)^2{\rm E}\left[x^2\right]_{p(x)}+\left(\frac{\partial f}{\partial y}\right)^2{\rm E}\left[y^2\right]_{p(y)}+2\mu\frac{\partial f}{\partial y}{\rm E}\left[y\right]_{p(y)}+2\mu\frac{\partial f}{\partial y}{\rm E}\left[y\right]_{p(y)}+2\frac{\partial f}{\partial x}\frac{\partial f}{\partial y}{\rm E}\left[xy\right]_{p(x)p(y)}\\ &=&\mu^2\cdot 1+\left(\frac{\partial f}{\partial x}\right)^2\cdot\frac{1}{4}+\left(\frac{\partial f}{\partial y}\right)^2\cdot\frac{1}{4}+2\mu\frac{\partial f}{\partial y}\cdot 0+2\mu\frac{\partial f}{\partial y}\cdot 0 + 2\frac{\partial f}{\partial x}\frac{\partial f}{\partial y}\cdot 0\\ &=&\mu^2+\frac{1}{4}\left(\left(\frac{\partial f}{\partial x}\right)^2+\left(\frac{\partial f}{\partial y}\right)^2\right)\tag{10} \end{eqnarray}$

以上で、深度値の ${\rm M}_1,{\rm M}_2$ が求まりました。(近似ですがね。)

機械学習基礎理論独習

誤りがあればご指摘いただけると幸いです。数式が整うまで少し時間かかります。リンクフリーです。

バリアンスシャドウマップのモーメント計算の行間を補う

はじめに

実際のアルゴリズム