Random Picking

特定の図形中のランダムな位置を取得したい。

棄却サンプリング

一番シンプルな方法は図形を全て含む領域で一様にサンプリングした後に、実際に図形内に含まれていれば採用、含まれていなければ棄却とする棄却サンプリングを使うことだ。しかし、

  • 図形の内外判定を行う必要がある。
  • 棄却された乱数の生成コストは丸々無駄になる。
  • 最終的に棄却される数は事前にわからないので生成する乱数の数を事前に決められない。

などいくつかのデメリットがある。棄却せずに直接マッピングできるならそれに越したことは無い。そのようなことが出来る図形についてまとめた。

平行四辺形

平行四辺形は平行な軸を二つ持つ図形なので、その二つの軸を基底ベクトル$\mathbf{e_x},\mathbf{e_y}$と見立て、独立な[0,1]の範囲の乱数$c_x,c_y$をその基底での座標とすることで内部の点を一様にサンプリングすることが出来る。

$ p = c_x \mathbf{e_x} + c_y \mathbf{e_y} $

MATLABでは以下のようになる。

{% highlight c++ %} n = 4096; cx = rand(n,1); cy = rand(n,1); ex = [0.5,0.0]; ey = [0.25,0.5]; x = cx * ex(1) + cy * ey(1); y = cx * ex(2) + cy * ey(2); scatter(x,y); {% endhighlight %}

parallelogram

三角形

平行四辺形は点対称に配置された二つの合同な三角形に分解できる。欲しい結果はそのうちの片方なので、まずは平行四辺形としてサンプリングし、仮に欲しいのとは逆側の三角形をサンプリングしてしまった場合($c_x+c_y > 1.0$)はそのサンプリング点の座標を平行四辺形の重心に対して$\pi$だけ回転すれば良い。重心$(c_x,c_y)=(0.5,0.5)$での点対称の座標は次のようになる。

$ \left( \begin{array}{ccc} c_x'\\\ c_y'\\\ \end{array} \right) = \left( \begin{array}{ccc} cos(\pi) & -sin(\pi)\\\ sin(\pi) & cos(\pi) \\\ \end{array} \right) \left( \begin{array}{ccc} c_x - 0.5\\\ c_y - 0.5\\\ \end{array} \right) + \left( \begin{array}{ccc} 0.5\\\ 0.5\\\ \end{array} \right) = \left( \begin{array}{ccc} 1.0 - c_x \\\ 1.0 - c_y \\\ \end{array} \right) $

MATLABでは以下のようになる。

{% highlight c++ %} n = 4096; cx = rand(n,1); cy = rand(n,1); e1 = [0.5,0.0]; e2 = [0.25,0.5];

for i = 1:4096 if cx(i)+cy(i) > 1.0 cx(i) = 1.0 - cx(i); cy(i) = 1.0 - cy(i); end end

x = cx * e1(1) + cy * e2(1); y = cx * e1(2) + cy * e2(2); scatter(x,y); {% endhighlight %}

triangle

正多角形

三角形が複数枚集まった形として正多角形がある。正多角形は合同な二等辺三角形から成るので次のような手順を踏むと正多角形内のサンプリングを行うことができる。

  1. 二等辺三角形の二つの等辺を基底とする座標系での座標を二つの乱数から決定する。
  2. 正N角形であれば、N個の三角形から成るのでどの三角形を選択するかをさらに一つの乱数から決定する。
  3. 決定された三角形の位置まで回転させる。

正六角形の場合をMATLABで行ったものが以下だ。

{% highlight c++ %} n = 4096; cx = rand(n,1); cy = rand(n,1); cz = rand(n,1); e1 = [3^(0.5)/2, 0.5]; e2 = [3^(0.5)/2, -0.5];

for i = 1:4096 if cx(i)+cy(i) > 1.0 cx(i) = 1.0 - cx(i); cy(i) = 1.0 - cy(i); end end

x = cx * e1(1) + cy * e2(1); y = cx * e1(2) + cy * e2(2); rot = double(mod(int8(cz * 6),6)) * pi / 3.0; x2 = cos(rot) .* x + -sin(rot) .* y; y2 = sin(rot) .* x + cos(rot) .* y; scatter(x2,y2); {% endhighlight %}

hex

円は、正多角形の三角の数を無限に増やした極限と見立てることができる。サンプリングの手順は次のようになる。

  1. $\theta$方向を向いている無限小の頂角を持つ二等辺三角形の二つの辺$\mathbf{e_x},\mathbf{e_y}$を基底ベクトルと見なして、三角形内のサンプリングで行ったのと同様に、それぞれの基底ベクトルの大きさを乱数値$c_x,c_y$としてサンプリング位置を決定する。
  2. 二等辺三角形の頂角は無限に小さいので、極座標系の動径$r$は単純に$c_x+c_y$となる。
  3. [0,1)の乱数$c_z$を$2\pi$でスケール、$[0,2\pi)$とし極座標の偏角$\theta$とする。

まとめると次のようになる。

$ r = \begin{cases} c_x + c_y & c_x + c_y \le 1 \\\ 1 - (c_x + c_y) & \text{otherwise} \end{cases} \\\ \theta = 2 \pi c_z $

これでも問題がないが、少し最適化を行うことができる。一様分布の確率分布(PDF)は次のようになっている。

$ f(t) = \begin{cases} 1 & 0 \le t < 1 \\\ 0 & \text{otherwise} \end{cases} \\\ $

ここでr、つまり二つの一様分布の和の確率分布関数(PDF)を考える。ある独立な二つの確率分布関数(PDF)があったとき、それらの和の確率分布は次のようになる。

$ f_{x+y}(t) = \int^{+\infty }_{-\infty } f_x(t-x) f_y(x) dx $

このことから二つ係数$c_x,c_y$の分布の和の分布は次のようになる

$ f_{c_x+c_y}(t) = \begin{cases} t & 0 \le t < 1 \\\ 2-t & \text{otherwise} \end{cases} \\\ $

累積分布関数(CDF)$F_{c_x+c_y}(t)$は確率分布関数(PDF)の積分なので結局次のようになる。

$ F_{c_x+c_y}(t) = \begin{cases} \cfrac{1}{2} t^2 & 0 \le t < 1 \\\ -\cfrac{1}{2}(t-2)^2 + 1 & \text{otherwise} \end{cases} \\\ $

あるCDFが与えられたときに、対応する確率変数はCDFの逆関数から算出されるので、結局最初の式のrは次のようになる。

$ r = c_x^{1/2} \\\ \theta = 2 \pi c_z $

$(r,\theta)$の算出に二つの乱数値($c_x,c_z$)でサンプリングができるようになった。

MATLABでは以下のようになる。

{% highlight c++ %} n = 4096; cx = rand(n,1); cy = rand(n,1); r = sqrt(cx); theta = cy * 2 * pi; x = cos(theta) .* r; y = sin(theta) .* r; scatter(x,y); {% endhighlight %}

circle

画像の濃淡を確率と見立てた場合

最後に少し毛色が違うが、濃淡画像で確率密度が表された場合のサンプリングについて書く。画像は内部的には一本の配列から構成されていることが多いので、その累積分布関数(CDF)を求めそれに対して、乱数値一つからサンプリングができる。しかしこれでは画素位置のみの疎なサンプリングになってしまう。

これの解決としては行毎に累積分布関数(CDF)を求め、さらに列用の累積分布関数(CDF)を一つ出す方法がある。

  1. 全ての行それぞれに対して、累積分布関数(CDF)を求める。
  2. 行毎に濃淡の合計値を出す。それをまとめた累積分布関数(CDF)を出し、これを列方向の累積分布関数(CDF)とする。
  3. サンプリングするときは、列の累積分布関数(CDF)から行を出す→最も近い行の累積分布関数(CDF)から列を出す。
  4. 対象の座標に対して何らかの補間処理(バイリニアなど)を施す。

詳しくは、Monte Carlo Rendering with Natural Illuminationの3.2を参照されよ。

その他メモ

  • 円への大雑把なマッピングとしてconcentric map(正式名称不明)がある。計算量は多少は少ないが円の内部で一様には分布しない。