この記事はレイトレアドベントカレンダー2020の初日の記事です。
モンテカルロ法では「乱数」を使って積分を数値的に行います。
その「乱数」は一体どのようなものを使えばいいのでしょうか？
早いしXorshiroとかでいいのでしょうか？
PRR bookを見てみると、なぜか丸ごと一章も使って解説しています。
ここはちゃんと頑張ると結果が良くなるところのようです。
頑張りたいところなのですが、指標が欲しいです。
レンダラーの「良さ」みたいなのを測るために、
レンダリング画像とリファレンス画像とのMSEをプロットするようなものもよく見ますが、
もっと「乱数」列それ自体に近づき観察して、それの「良さ」みたいなのを評価したくなります。
ここからは乱数と呼ばずにサンプル列と呼びます。
サンプル列は複数のサンプルから成り、サンプル一つ一つは何次元かの値になっています。
サンプル1 = (x0,y0,z0,w0,...)
サンプル2 = (x1,y1,z1,w1,...)
サンプル3 = (x2,y2,z2,w2,...)
...
違うサンプルの同じ次元の数値は、同じ積分変数を指しています。
評価方法として、まず思い浮かぶのは隣接する次元同士で「一様」に分布することを確かめることです。
どういうことかパストレで例えてみます。
「薄レンズモデルのカメラから見た、矩形のエリアライトの直接光に照らされる物質」
というシーンを考えてみます。
このときサンプルはどのように使われるかというと...
となりサンプルは合計7次元必要になります。
6次元目と7次元目は矩形のライトの上の位置のサンプルになるのですが、
例えばあるサンプル列で6次元目と7次元目が、次のようにプロットされたとします。
import matplotlib.pyplot as plt
import random
plt.axes().set_aspect('equal', 'box')
plt.scatter(
[random.random() for i in range(256)],
[random.random() for j in range(256)],
c=[range(256)], s=50)
plt.show()
この[0,1)^2の範囲をそのまま矩形のライト上の位置だと思ってください。
ぐっと睨むと、ちょっとこれだとまずそうだという気分になってきます。
例えば、点が全くない空白地帯や、逆に点が沢山重なっている場所に、
特徴的なもの(物体のエッジとか、テクセル境界とか)
が重なってしまうとその特徴をスキップしたり、余計に考慮してしまって収束が遠のいてしまいそうです。
Pythonのrandomを使っているので数学的には一様ではあるはずなのですが、
モンテカルロ法をしているときに欲しい「一様」とはちょっと違いそうです。
代わりに、もしこんな感じにサンプルされていたらどうでしょうか。
import math
import matplotlib.pyplot as plt
plt.axes().set_aspect('equal', 'box')
plt.scatter([(0.5 + 0.628706721037808 * i) % 1.0 for i in range(256)],
[(0.5 + 0.733891856627126 * i) % 1.0 for i in range(256)],
c=[range(256)], s=50)
plt.show()
とてもよさそうです！
このようにプロットすることで、欲しい「一様さ」が得られているか
(そもそも実装ミスをしていないかなども)を簡単に確認することができます。
ではできた二次元のサンプル列を、コピーして他の次元にもコピーして使いまわしましょう！めでたしめでたし！
....とは残念ながらなりません。
よくできた二次元のサンプル列があるならばそれを異なる次元で使いまわした方が、処理的に有利そうです。
しかしそれはやってはいけません。なぜならば
モンテカルロ積分は2次元の積分が複数集まったもの $f(x0,x1,x0,x1,...)$
ではなく、N次元の一つの積分 $f(x0,x1,x2,x3...)$だからです。
例え隣り合っていない次元だとしても、それは互いに影響を与え合っては(相関しては)いけません。
先のパストレの例を再掲します。
何らかの使いまわしの結果、2次元目と7次元目が同じ値になるとします。
この場合はサブピクセルとエリアライトのサンプルに相当します。
「サブピクセルの右側の領域を使うと必ず、エリアライトも右半分がサンプルされる」
といったことが起こってしまい、これは紛れもなくバイアスです。
この話は、全く同じ値を使いまわさなくても起こる事に注意してください。
「i次元目をちょっといじってj次元目を作った結果、(i,j)次元間に相関が起きた」
場合も同じようなバイアスが乗ってしまいます。この辺りの話は[1]にもまとめられています。
ということで、N次元のサンプルをプロットする際は
(1,2),(2,3),(3,4),...(N-1,N) のN-1個のプロットではなく、
N個の次元から2個の組み合わせ $ {}_n \mathrm{ C }_2 = N(N-1)/2 $ のプロットが必要だとわかりました。
実際にちゃんと次元の全組み合わせをプロットしないとダメな場合を例示します。
以下のコードはあるサンプラーを実装し、プロットしています。
import math
import matplotlib.pyplot as plt
def ri(n, d):
base = [2,3,5,7,11,13,17,19,23,29][d]
nume = 0
deno = 1
while n > 0:
nume = nume * base + (n % base)
deno = deno * base
n //= base
return float(nume) / float(deno)
fig = plt.figure(figsize=(8.0, 10.0))
for d0 in range(4):
for d1 in range(d0+1,5):
ax = fig.add_subplot(5, 4, d0 + d1 * 4 + 1)
ax.set_aspect('equal', 'box')
ax.tick_params(
labelleft=False, labelbottom=False,
bottom=False, left=False,right=False,top=False)
ax.scatter([ri(i,d0) for i in range(128)],
[ri(i,d1) for i in range(128)],
c=[range(128)], s=10)
if d0 == 0:
ax.set_ylabel(str(d1+1))
if d1 == 4:
ax.set_xlabel(str(d0+1))
plt.show()
縦と横に振られた数字が次元を表しています。
例えば、一番左上のプロットは、「1次元と2次元」のプロットになります。
先ほど注目していたのは「隣接する次元」だけなので、この図全体でいうところの、
対角線上にある(1,2),(2,3),(3,4),(4,5)次元のプロットになります。
対角線上のプロットをみても特に問題が無さそうに見えます。
しかし左下の(1,5)次元などは何やら怪しい模様が出てきており、
積分した結果にも何らかのアーティファクトが出ることが予想されます。
全ての次元の組み合わせをプロットしたほうがよさそうです。
全ての次元の組み合わせをプロットするとして、あとは目で「良さ」を評価すればいいのでしょうか？
ぱっと見の良さが肉薄している次の二つのサンプルの「良さ」を評価してくださいと言われたら困ってしまいそうです。
import math
import matplotlib.pyplot as plt
fig = plt.figure(figsize=(8.0, 10.0))
ax0 = fig.add_subplot(1, 2, 1)
ax1 = fig.add_subplot(1, 2, 2)
ax0.set_aspect('equal', 'box')
ax1.set_aspect('equal', 'box')
ax0.scatter([(0.5 + 0.628706721037808 * i) % 1.0 for i in range(256)],
[(0.5 + 0.733891856627126 * i) % 1.0 for i in range(256)],
c=[range(256)], s=50)
ax1.scatter([ri(i,2) for i in range(256)],
[ri(i,3) for i in range(256)],
c=[range(256)], s=50)
plt.show()
答えから書くと、サンプル列のパワースペクトルを見て判断されることが多いようです。
[2]を参考に実装します。
import math
import numpy as np
from PIL import Image, ImageChops
from tqdm import trange
# Spectrum表示の画像サイズ
# 本当はサンプル数に合わせてサイズを変えなければいけないが、
# 利便性を取ってユーザー設定できるように。
SPECTRUM_IMAGE_SIZE = 256
#
def computePowerSpectrum2d(samples, f, d0, d1):
ds = 2.0 * math.pi * (f[0] * samples[d0] + f[1] * samples[d1])
st = np.sum(np.sin(ds))
ct = np.sum(np.cos(ds))
return (st * st + ct * ct) / len(samples[0])
#
def drawPowerSpectrum(samples, d0, d1):
#img = Image.new('RGB', (SPECTRUM_IMAGE_SIZE,SPECTRUM_IMAGE_SIZE))
#pixels = img.load()
pixels = []
for y in range(SPECTRUM_IMAGE_SIZE):
line = []
for x in range(SPECTRUM_IMAGE_SIZE):
f = [
(x/SPECTRUM_IMAGE_SIZE-0.5) * 2.0 * SPECTRUM_IMAGE_SIZE,
(y/SPECTRUM_IMAGE_SIZE-0.5) * 2.0 * SPECTRUM_IMAGE_SIZE
]
pv = computePowerSpectrum2d(samples, f, d0, d1)
# tone mapping
pvi = int(math.log2(pv + 0.25) * 255)
line.append(pvi)
#pixels[x,y] = (pvi, pvi, pvi)
pixels.append(line)
return pixels
余談ですが、当初はPythonベタ書きだったのですが、ここの処理が遅すぎて、
途中でNumPyに書き換えたところ、途端に100倍近く高速になって驚きました。
コア数分位は早くなってくれるだろうと思ったのですが、それ以上でした。
実装できたものを使って、試しに先ほどの二つのサンプル列に対して適用してみます。
import math
import matplotlib.pyplot as plt
img0 = drawPowerSpectrum(
np.array(
[[(0.5 + 0.628706721037808 * i) % 1.0 for i in range(2048)],
[(0.5 + 0.733891856627126 * i) % 1.0 for i in range(2048)]]),0,1)
img1 = drawPowerSpectrum(
np.array(
[[ri(i,2) for i in range(2048)],
[ri(i,3) for i in range(2048)]]),0,1)
fig = plt.figure(figsize=(8.0, 10.0))
ax0 = fig.add_subplot(1, 2, 1)
ax1 = fig.add_subplot(1, 2, 2)
ax0.set_aspect('equal', 'box')
ax1.set_aspect('equal', 'box')
ax0.pcolorfast(img0, cmap="plasma")
ax1.pcolorfast(img1, cmap="plasma")
plt.show()
ぱっと見では同じように見えた、サンプル列も、全く違う周波数特性を持つことがわかりました。
ここまでのまとめです。
最後にこれまでの実装を一つにまとめましょう。
[3]の表示の仕方をパクります。
import matplotlib.pyplot as plt
import numpy as np
from tqdm.notebook import trange, tqdm
# プロット時の最大のサンプル数
MAX_PLOT_SAMPLE = 128
def getSamples():
return np.array(
[
[ri(i,0) for i in range(2048)],
[ri(i,1) for i in range(2048)],
[ri(i,2) for i in range(2048)],
[ri(i,3) for i in range(2048)],
[ri(i,4) for i in range(2048)],
[ri(i,5) for i in range(2048)],
])
#
def plotStairs(samples):
numDim = len(samples)
numSample = len(samples[0])
dims = []
for d0 in range(numDim-1):
for d1 in range(d0+1, numDim):
dims.append([d0,d1])
fig = plt.figure(figsize=(12.0 * (numDim+1)/numDim, 12.0))
plt.tight_layout() # きいていない？
for dim in tqdm(dims):
d0 = dim[0]
d1 = dim[1]
axPlot = fig.add_subplot(numDim, numDim+1, d0 + 1 + (d1-1) * (numDim+1))
# [0,1]に限定する
axPlot.set_xlim(0,1)
axPlot.set_ylim(0,1)
#
major_ticks = np.arange(0, 1, 1/4)
axPlot.set_xticks(major_ticks)
axPlot.set_yticks(major_ticks)
axPlot.grid(which='both')
axPlot.tick_params(
labelleft=False, labelbottom=False,
bottom=False, left=False,right=False,top=False)
#
axPlot.set_aspect('equal', 'box')
# 縦と横にラベル
if d0 == 0:
axPlot.set_ylabel(str(d1+1))
if d1 == numDim-1:
axPlot.set_xlabel(str(d0+1))
axPlot.scatter(samples[d0][:MAX_PLOT_SAMPLE], samples[d1][:MAX_PLOT_SAMPLE], c=range(MAX_PLOT_SAMPLE), s=5)
#
axPower = fig.add_subplot(numDim, numDim+1, d1 + 1 + d0 * (numDim+1))
axPower.tick_params(
labelleft=False, labelbottom=False,
bottom=False, left=False,right=False,top=False)
if d0 == 0:
axPower.xaxis.set_label_position("top")
axPower.set_xlabel(str(d1+1))
if d1 == numDim-1:
axPower.yaxis.set_label_position("right")
axPower.set_ylabel(str(d0+1))
img = drawPowerSpectrum(samples, d0, d1)
axPower.pcolorfast(img, cmap="plasma")
plt.show()
#
plotStairs(getSamples())
ついにできました！
このツールを駆使して、君だけの最強のサンプラーを作ろう！
Jamorn Sriwasansakさんに、実装について教わりました。
[1] E. Turquin From Ray to Path Tracing: Navigating through Dimensions
[2] T. Schlömer, O. Deussen Accurate Spectral Analysis of Two-Dimensional Point Sets Journal of Graphics, GPU, and Game Tools, 15(3):152−160, 2011
[3] W. Jarosz, A. Enayet, A. Kensler, C. Kilpatrick, P. Christensen. Orthogonal array sampling for Monte Carlo rendering. Computer Graphics Forum (Proceedings of EGSR), 38(4):135–147, July 2019.
ここでは使いませんでしたが、実際にツールとして使う時の便利関数です。
CSVファイルに記述されたサンプル列をロードします。
import numpy as np
# samples[dim][id]の形になっていることに注意
def readSamples(filePath):
# 最大の次元
MAX_DIMENSION = 5
# 最大のサンプル数
MAX_SAMPLE = 2048
#
with open(filePath) as f:
header = f.readline().split()
numSample = min(int(header[0]), MAX_SAMPLE)
numDim = min(int(header[1]), MAX_DIMENSION)
samples = np.empty([numDim, numSample])
for sn in range(numSample):
sampleStr = f.readline().split()
for di in range(numDim):
samples[di][sn] = float(sampleStr[di])
return samples