サンプリングとレンダリングとガチャについて

# はじめに
できるだけ間違いがないように努力はしますが, 私は数学が本当に苦手なため, 必ず自分で考えるようにしてください.

前半は, 3D CGの積分におけるサンプリングについて調べたこと, 後半は, ゲームにおけるガチャについて調べたことになります.

3D CGのレンダリングは積分問題だということができます. これは, レンダリング方程式を見ると明らかです.

$L_{o}(x\omega_{o},\lambda,t)=L_{e}(x,\omega_{o},\lambda,t) + \int_{\Omega}f_{r}(x,\omega_{i},\omega_{o}\lambda,t)L_{i}(x,\omega_{i},\lambda,t)(\omega_{i}\cdot n)d\omega_{i}$

レンダリングの積分は普通, 解析的に解けずモンテカルロ法等を使って解きます.
モンテカルロ法は乱数を用いて統計的に数値計算を行う手法です. ここでは, 乱数を用いて統計的にレンダリング結果を推定することです.

基本的なことについて書いておきます.

## 疑似乱数
レンダリングにおけるモンテカルロ法やゲームのガチャでは, 速度的理由から擬似乱数を使用します. 一般的に乱数は局所的な分布について何も保障していません. つまり, 数千や数万のオーダーでは一様に分布しているが, 数十や数百のオーダーでは偏りがあるかもしれません.

## モンテカルロ法
モンテカルロ法は, 乱数を用いて標本抽出をする方法と言い換えることができます. 完全な乱数の場合, サンプル数$N$にたいして$N^{-1}$の速度で解に収束するそうです.

## 準モンテカルロ法
乱数のようにみえるが, 規則的な構造と順序を持っている点列を用いる方法です. この点列は局所的にも偏りがありません.
問題とサンプル点の生成方法によって劇的に解への収束速度が速くなります. 常にモンテカルロ法より収束が速くなるとは限りません. 相関のないサンプルに使用すべきではありません.

## ブルーノイズ
正確な定義と少し違うらしいのですが, CG業界ではホワイトノイズと異なり, いいかんじに散らばった点群のこと, ポアソンディスクサンプリングのように一定以上離れた点群をいいます[(Wikipedia カラードノイズ)](https://ja.wikipedia.org/wiki/%E3%82%AB%E3%83%A9%E3%83%BC%E3%83%89%E3%83%8E%E3%82%A4%E3%82%BA).
他に良い記事があるためここでは扱いません. "Capacity-Constrained Voronoi Tessellation"を疑似コードそのままに実装したのですが, かなり計算に時間がかかります.

## 逆関数法
ある確率分布に従うランダムな値を得たいことがあります. 3D CGの文脈では重点的サンプリングが例になります.

ある確率分布$P(x)$に従う乱数を, 一様分布に従う乱数から生成する方法です.
$U$が一様分布に従うとすると,
$$P(U \le u) = u \  (0 \le u \le 1)$$
$P(x)$の累積分布関数を$F(x)$とすると,
$$U \le u \Leftrightarrow F^{-1}(U) \le F^{-1}(u)$$
が成り立つため, 次も成り立つ.
$$P(F^{-1}(U) \le F^{-1}(u)) = u$$
$F^{-1}(u)=x$とおくと,
$$P(F^{-1}(U) \le x) = F(x)$$
$F^{-1}(U)$は累積分布関数が$F(x)$となる確率分布$P(x)$に従う.

積分が解析的に解ける場合, 上記の通りに計算すればよいのですが, 現実の問題はそう簡単ではありません. とは言え効率を考えなければ頭を使わず解くことができます.

ものすごく単純なゲームのガチャ, 3枚のカードが[0.1, 0.4, 0.5]の確率で出るガチャを考えます. 確率を先頭から累積したテーブル[0.1, 0.5, 1.0]を作り, サイコロを振ります. このように正規化している場合 [0 1)の値が出るサイコロを振ります. サイコロの結果が越える場所を探索すればよいでしょう. サイコロの結果が0.5なら2番目のカードが選ばれます.

# CGにおけるサンプリング手法
モンテカルロ法で積分を解くことになりますが, より収束が速くなるようにさまざまな手法が研究されています.
次に, 調査したものを羅列していきます. 私はエンジニアであって, 研究者ではないので, 網羅性はありません. Siggraph 2019での"My favorite Samples"([[1]](#[1]))から始めてみるのがいいと思います. 
ピクセルをスーパーサンプリングする場合を想定して, 2次元に1024点プロットした図を添付しています.

##手法
### 乱数
ホワイトノイズはこうなります. どの周波数も均等に出現するはずなので, 小さな空地から大きな空地まで均等に出現します.
![random_points.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/27065/a529ae45-6041-a5ff-bde5-2abf63eeda28.png)

### 層化抽出法(stratified sampling)
一様格子で分割し, 格子内でランダムにサンプルします. 次元数のべき乗のサンプル数が事前にわかっている方がいいです. 格子（でなくてもいいですが）の面積を変えて簡易的に密度を変えることができます.

![stratified_points.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/27065/084e1376-4073-bbd0-9c36-26b78c34a5b6.png)

### Halton列
2次元では4象限を順繰りに出現する点列を作成できます([[2]](#[2])).
1024点程度では知覚できないのですが, もっと多くすると構造を持っていることがわかります.

![halton_points.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/27065/5c71978c-bf03-6f4b-8516-cca1a24009ba.png)

### Hammersley点群
はじめにサンプル数を決める必要があります.
![hammersley_points.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/27065/333d9ea9-865c-5271-01b8-ceef5afd2d4c.png)

### Rank1 Lattice
パラメータの決定が難しいですが, 2次元だと総当たりでもそれほどの時間はかかりません([[4]](#[4]), [[5]](#[5])). 一様格子と同じく非常に整った構造のため, 意図的に苦手なパターンを作ることは簡単だと思います.
![rank1lattice_points.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/27065/c9c0e5ee-9a0e-f2a4-bc76-018d9aac04e0.png)

### Progressive Jittered Sample Sequences
4の倍数でないとならないという制約はありますが, progressiveにサンプル点を生成できます([[6]](#[6])).
Multi-Jitteredもありますが, 次に点を置くべき空き領域探索を単純な乱択で行うため, 計算量が爆発します. 2分木で空き領域探索を高速化する手法が提案されています([[7]](#[7])).
![progressive_jittered_points.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/27065/1571834f-24cb-5f08-12a4-2e8183e85490.png)

### R2
[[6]](#[6])を基に, 黄金比を使って簡易化したものということでいいのかな？([[8]](#[8]))
R2, Jittered R2, Random Jittered R2のバリエーションです.

![r2_points.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/27065/c64be43e-2083-b529-58c9-63a269639c50.png)

![jittered_r2_points.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/27065/7b295031-810c-37a3-581b-875735301986.png)

![random_jittered_r2_points.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/27065/d521402d-8f3e-67f8-e120-0b9411649953.png)

## 実験
$[-0.5\ 0.5]\times[-0.5\ 0.5]$で$e^{-(x^2 + y^2)}$のモンテカルロ積分を行いました([octave-online](https://octave-online.net/)).

```
[x,y] = meshgrid(-0.5:.05:0.5);                                
z = exp(-(x.^2 + y.^2));
surf(x,y,z)
```

およそ$0.85112$になるはずです.
次のグラフは, サンプル数を横軸に, 二乗平均平方根誤差rmseを縦軸にしています.
横軸は常用対数です. Jittered R2がないのは, 計算が終わらないからです.

![integration00.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/27065/2fd64213-0ea8-658b-d93f-f9d2cf41c3b5.png)

右にいくほどサンプル数が増加し, エラーが減少しています. 乱数は一番上のラインですから, どの手法も乱数より収束が速いといえます.
連続な関数$e^{-(x^2 + y^2)}$の積分をした場合ですので, 積分対象によって結果は変わります.

## 結果
乱数は局所的な分布については何も保障していないよ, ということがわかっていただければいいと思います.
どれが最高とか言えないのは, 適材適所でしかないからです.
乱数やStratified以外は, 構造や規則があるため, その知識を使って上手いフィルタリングができたりするのでしょうか.

# ガチャの実装について
前半と全く異なる話題ですが, いろいろ事情があるのです. ほぼポエムなので飛ばしてください.

実装次第でガチャは局所的に偏ります. レアリティを先ず選択し, レアリティ内でカードを選択する2段式だと偏ると思います. "ガチャ 実装"などのキーワードで調べると, 2つぐらいこのような実装を見つけました.

問題は1段目のレアリティを選択する部分です. コンピュータで疑似乱数で少数の母集団から乱択するとき, 局所的に偏ります. レアリティは高々5個くらいでしょうから, 時間的に偏ると思います.
ほぼ間違いなく, サーバは複数のプレーヤでサイコロを共有する実装でしょうから, 1回ガチャは偏りが起きないと思います. 先の2段式だと, 最高レアリティが出るときはやたら出たり, 出ないときはでなかったり. 10連では, 出るときは2, 3枚同時に出ることが頻発するということです.
実験してみたところ, 10連を10万回繰り返して偏りあるなぁ, なので気にするほどでもないのですが.

## 対策？
普通に律儀に逆関数法をするといいと思います. あなたの考えた素晴らしいアルゴリズムは, おそらく偏っています.
やりすぎとは思いますが, 確率を累積したテーブルを作成する前に, 毎回シャッフルするといいです.

# まとめ
乱数（ホワイトノイズ）の局所的な分布は, あなたの感覚と違うかもしれませんよ, ということだけ覚えてもらえたらいいです.

# 参照
<a name="[1]">[1]</a> "My favorite Samples", ACM Siggraph 2019 Course,
https://sites.google.com/view/myfavoritesamples
<a name="[2]">[2]</a> John H Halton, "Algorithm 247: Radical-inverse quasi-random point sequence"
<a name="[3]">[3]</a> Hammersley Points on the Hemisphere,
http://holger.dammertz.org/stuff/notes_HammersleyOnHemisphere.html
<a name="[4]">[4]</a> Sabrina Dammertz, Alexander Keller,
"Image Synthesis by Rank-1 Lattices",
https://link.springer.com/chapter/10.1007/978-3-540-74496-2_12
<a name="[5]">[5]</a>
Program on Quasi-Monte Carlo and High-Dimensional Sampling Methods for Applied Mathematics Opening Workshop, Lattice Rules for Quasi-Monte Carlos - Pierre L’Ecuyer, Aug 28, 2017
https://www.slideshare.net/SAMSI_Info/program-on-quasimonte-carlo-and-highdimensional-sampling-methods-for-applied-mathematics-opening-workshop-lattice-rules-for-quasimonte-carlos-pierre-lecuyer-aug-28-2017
<a name="[6]">[6]</a> Per Christensen, Andrew Kensler, Charlie Kilpatrick,
"Progressive Multi-Jittered Sample Sequences",
EGSR 2018,
https://graphics.pixar.com/library/ProgressiveMultiJitteredSampling/paper.pdf
<a name="[7]">[7]</a> Matt Pharr,
"Efficient Generation of Points that Satisfy Two-Dimensional Elementary Intervals",
JCGT 2019,
http://jcgt.org/published/0008/01/04/
<a name="[8]">[8]</a>
http://extremelearning.com.au/a-simple-method-to-construct-isotropic-quasirandom-blue-noise-point-sequences/
<a name="[9]">[9]</a> Colas Schretter, Leif Kobbelt, Paul-Olivier Dehaye,
"Golden Ratio Sequences for Low-Discrepancy Sampling",
JGT 2012,
https://www.graphics.rwth-aachen.de/publication/2/jgt.pdf