モーメントシャドウマップについて
はじめに
分散シャドウマップに続いてモーメントシャドウマップについてです。
モーメントシャドウマップの方が少し難しい感じでなかなか理解がすすまなかったんですが、論文のほかに補足的なドキュメントが出ていてそちらを見ていたら大体の流れが理解できてきた気がするのでまとめようと思います。
アルゴリズムの部分で出てきた式をなんとなく追いかける目的で書いているため色々と細かい部分は飛ばしています。
参考にしたもの
研究者の方のページに論文やスライド、発表動画、補足ドキュメントなど情報が色々載っていました。
導出について詳しく書いてあった補足のドキュメントはこの pdf です。
http://cg.cs.uni-bonn.de/aigaion2root/attachments/MomentShadowMappingSupplementary.pdf
モーメントシャドウマップ
まずモーメントシャドウマップについて説明してみます。
分散シャドウマップと少し似ていて改良した方法になっています。
分散シャドウマップでは二次のモーメントまでを使いますが、モーメントシャドウマップではより高次のモーメント (通常は四次) まで利用します。
ライトブリーディングという遮蔽物の境界で本来影であるはずのところに光の筋ができてしまう問題が改善されているようです。
MIPMAP やガウシアンブラーなどのアンチエイリアシング手法も同様に使用できます。
アルゴリズム
論文の Algorithm 2 にニの倍数のモーメントでの計算方法が書かれています。
先に使われている変数をいくつか紹介します。
$z$ はライトからの深度です。
$\textbf{b}$ は深度からモーメントの列ベクトルを返す関数です。
$$
\textbf{b}(z) := (z, z^{2}, \ldots, z^{n})^{T}
$$
$b$ は $\textbf{b}(z)$ の平均値です。
$\hat{\textbf b}$ は0乗からn-1乗までのモーメントの列ベクトルを返す関数です。
$$
\hat{\textbf b}(z) := (1, z, z^{2}, ..., z^{n-1})^{T}
$$
$\hat{b}$ は $\hat{\textbf b}(z)$ の平均値です。
Algorithm 2
- $$n := \frac{m}{2} + 1$$
モーメントの次数$m$から求める深度の数$n$が決まります。
- $$ B := (b _ {j+k-2}) _ {j,k=1} ^ n $$
モーメント平均値から$B$行列を作ります。
- $$ z _ 1 := z _ f $$
深度のうち一つは描画する部分の深度とします。
- $$ c := B^{-1} \cdot \hat{\textbf b}(z_{1})$$
$B$の逆行列とベクトル $\hat{\textbf b}(z_{1})$ の積を計算します。これを$c$とします。
- $$ \sum _ {k=1} ^ n c _ k z ^ {k-1} = 0 $$
$c$と深度の $n-1$ 次の方程式を解いて $n-1$ 個の深度を求めます。
- $$\hat{A} := (z _ i ^ {j-1}) _ {j,i=1} ^ n $$
求めた深度を使って $\hat{A}$ 行列を求めます。この行列はヴァンデルモンドマトリックスと呼ばれ逆行列を作れます。
- $$ w = \hat{A}^{-1} \hat{b} $$
$\hat{A}$ の逆行列とモーメント平均値のベクトルから各深度の比率を求めます。
- $$ G := \sum_{i=1,z _ i \lt z _ f} ^ n w _ i $$
描画部分の深度よりも小さい深度の比率を足して影の濃さを計算します。
導出
モーメントシャドウマップでは$n$個の深度に集中して分布するとして深度と比率を求めて解いています。
モーメント平均値は次のように表せます。$w _ i$ は深度 $z _ i$ の分布比率です。
$$
\hat{b} = \sum _ {i=1} ^ n w _ i \hat{\textbf b}(z _ i)
$$
次に $\hat{\textbf b}(x _ {i})$ の直積の平均値を考えます。 これは$B$と等しくなります。
$$
\begin{aligned}
\sum _ {i=1} ^ n w _ i \hat{\textbf b}(z _ i)\otimes\hat{\textbf b}(z _ i) &= (\sum _ {i=1} ^ n w _ i z _ i ^ {j-1} z _ i ^ {k-1}) _ {j,k=1} ^ n \\
&= (\sum _ {i=1} ^ n w _ i z _ i ^ {j+k-2}) _ {j,k=1} ^ n \\
&= (b _ {j+k-2}) _ {j,k=1} ^ n \\
&= B
\end{aligned}
$$
変形して次の形にできます。
$$
\begin{aligned}
B &= (\sum _ {i=1} ^ n w _ i z _ i ^ {j-1} z _ i ^ {k-1}) _ {j,k=1} ^ {n} \\
&= (\sum _ {i=1} ^ {n}\sum _ {l=1} ^ {n} w _ {i} \delta _ {il} z _ {i} ^ {j-1} z _ {l} ^ {k-1}) _ {j,k=1} ^ {n} \\
&= \hat{A}\text{diag}(w)\hat{A} ^ {\textrm{T}} \\
&=
\begin{pmatrix}
1 & 1 & \cdots & 1 \\
z _ {1} & z _ {2} & \cdots & z _ {n} \\
\vdots & \vdots & & \vdots \\
z _ {1} ^ {n-1} & z _ {2} ^ {n-1} & \cdots & z _ {n} ^ {n-1}
\end{pmatrix}
\begin{pmatrix}
w _ {1} & & & 0 \\
& w _ {2} & & \\
& & \ddots & \\
0 & & & w _ {n}
\end{pmatrix}
\begin{pmatrix}
1 & z _ {1} & \cdots & z _ {1} ^ {n-1} \\
1 & z _ {2} & \cdots & z _ {1} ^ {n-1} \\
\vdots & \vdots & & \vdots \\
1 & z _ {n} & \cdots & z _ {n} ^ {n-1}
\end{pmatrix} \\
\hat{A} ^ {-1} B \hat{A} ^ {-\textrm{T}} &= \text{diag}(w)
\end{aligned}
$$
$\hat{A} ^ {T} B ^ {-1} \hat{A}$ も対角行列となります。
$$
\hat{\textbf b}(z _ {j}) ^ {\textrm{T}} B ^ {-1} \hat{\textbf b}(z _ {k}) =
\begin{cases}
0 & \text{if} & j \ne k \\
1 / w _ {j} & \text{if} & j = k
\end{cases}
$$
ここから Algorithm 2 のステップ4とステップ5の式が出てきます。
ステップ7の各深度の分布比率はモーメント平均値の式から出てきます。
$$
\begin{aligned}
\hat{b} &= \sum _ {i=1} ^ {n} w _ {i} \hat{\bf{b}}(z _ {i}) \\
&= \hat{A} w \\
w &= \hat{A} ^ {-1} \hat{b}
\end{aligned}
$$
$m=2$ の場合
$m=2$ の場合は分散シャドウマップと同じ式になります。
(分散シャドウマップの論文に出ていたのは影にならない確率なので実際には1から引いた値と同じ)
$$
\begin{aligned}
m &= 2 \\
n &= 2 \\
B &= \begin{pmatrix}
1 & b _ {1} \\
b _ {1} & b _ {2}
\end{pmatrix} \\
A &= \begin{pmatrix}
1 & 1 \\
z _ {1} & z _ {2}
\end{pmatrix} \\
z &= \frac{b _ {1} z _ {f} - b _ {2}}{z _ {f} - b _ {1}} \\
G &= \frac{(z _ f - \mu) ^ 2}{\sigma ^ 2 + (z _ f - \mu) ^ 2}
\end{aligned}
$$
$m=4$ の場合
Algorithm 3 の Hamburger 4MSM アルゴリズムというのが多分実際に実装に使われるアルゴリズムだと思います。
$B$ の逆行列が必要ですが計算誤差で正則でなくなってしまうことを避けるため$b$にバイアスを加えた $b'$ を使っているようです。
バイアスに使われている $\alpha$ は$0$以上の十分小さな値 (例えば $3.0 \cdot 10 ^ {-5}$ ) です。
$$
\begin{aligned}
m &= 4 \\
n &= 3 \\
b' &= (1 - \alpha) b + \alpha \begin{pmatrix}0.5 \\ 0.5 \\ 0.5 \\ 0.5\end{pmatrix}\\
B &= \begin{pmatrix}
1 & b' _ 1 & b' _ 2 \\
b' _ 1 & b' _ 2 & b' _ 3 \\
b' _ 2 & b' _ 3 & b' _ 4
\end{pmatrix} \\
A &= \begin{pmatrix}
1 & 1 & 1 \\
z _ {1} & z _ {2} & z _ {3} \\
z _ {1} ^ 2 & z _ {2} ^ 2 & z _ {3} ^ 2
\end{pmatrix} \\
G &= \begin{cases}
0 & \text{if} & z _ f \leq z _ 2 \leq z _ 3 \\
\frac{z _ f z _ 3 - b' _ 1 (z _ f + z _ 3) + b' _ 2}{(z _ 3 - z _ 2)(z _ f - z _ 2)} & \text{if} & z _ 2 < z _ f \leq z _ 3 \\
1 - \frac{z _ 2 z _ 3 - b' _ 1 (z _ 2 + z _ 3) + b' _ 2}{(z _ f - z _ 2)(z _ f - z _ 3)} & \text{if} & z _ 2 \leq z _ 3 < z _ f \\
\end{cases}
\end{aligned}
$$
$G$ は Wolfram alpha で同じ感じの答えが出るのが計算できました。
Simplify[Inverse[{{1,1,1},{z_f,z_2,z_3},{z_f^2,z_2^2,z_3^2}}][[2]].{1,b_1,b_2}]
Simplify[Inverse[{{1,1,1},{z_f,z_2,z_3},{z_f^2,z_2^2,z_3^2}}][[2;;3]].{1,b_1,b_2}.{1,1}]
思ったこと
最初は数か所の深度だけに分布すると仮定して計算するのは、それでいいのかなという疑問を持ちました。
ただ $m=2$ で分散シャドウマップで使われているチェビシェフの不等式と同じ結果となったのと、モーメントのみから正確な深度の分布を求めることはできないので下限の確率を求めるということがこの計算の意味なのかなと考えるとなんとなく納得できました。
ライトブリーディングが起こるような二箇所の深度に集中して分布する場合などにこの方法で対応できるということも理解できます。
ちょっと複雑だなという気はやっぱりしてしまいます。ライトブリーディングの問題が解決できれば他の方法でもいいかもとも思いました。
数式の表示
MathJax を使うようにしてみました。
同じ数式の書き方で VS Code などで Markdown 表示の確認ができて、はてなブログの Web 上で確認するよりも大分すばやく確認できて書きやすいと思います。
しかしはてなブログ特有の Markdown 記法と干渉したりして結局そのまま使えるわけではないので微妙のような気もします。
別のページに何か影響でてしまっているかとかもあまり確認してないので後で色々調整必要かもしれません。
