二次元の確率密度分布$f(x,y)$が与えられているとする．三次元空間$(x,y,z)$内に一様分布の乱数点$(x_i,y_i,z_i)$を生成し，次式を満たす箇所のみ点$(x_i,y_i)$を残すことで目的の乱数が得られる．一様分布を生成する$z$の範囲は$f(x,y)$の最大値$z_{\rm max}$を用いて$0 \leq z < z_{\rm max}$とすれば，$0\leq z < 1$としたときと等価な乱数が得られる．
\begin{eqnarray} z_i < f(x_i,y_i) \end{eqnarray}準備
必要なパッケージのimportとよく使う関数を定義しておく． plt.imshowとplt.savefigで画像を保存しても良かったけど画像サイズが変わってしまうのでcv2.imwriteを使うことにした．
import cv2
import numpy as np
def jet(img):
img_norm = np.zeros_like(img)
cv2.normalize(img, img_norm, 0, 255, cv2.NORM_MINMAX);
img_jet = cv2.applyColorMap(img_norm.astype(np.uint8), cv2.COLORMAP_JET)
return img_jet
確率密度分布を画像で表す
まずは，二次元の確率密度分布$f(x,y)$を画像として表す． 単純な発想では顕著性マップは注視位置の確率分布と考えられるので，適当な画像から顕著性マップ作ってみる． ここでは結果だけ載せておくけど使用したのはDeepGaze II．入力画像はDeepGaze IIの入力画像サイズに合わせて1024$\times$576を1024$\times$768にしてあるので少し縦に伸びている． ここでは入力画像も顕著性マップも何でもいいのでこれで良しとする．
sal = np.load("keiba_sal.npy")
norm = np.linalg.norm(sal, axis=None)
sal = sal/norm
cv2.imwrite("keiba_sal.png", jet(sal))
一様乱数の生成
$(x,y,z)$の各変数で独立に一様乱数を発生させれば$(x,y,z)$の一様乱数ができるけど，$(x,y)$に関しては画像の座標なので整数値に限定した． ここでは10000個の乱数点を発生させる．といっても，x,y,zにそれぞれ座標（乱数列）が代入されるだけ．
max = np.max(sal) n = 10000 x = np.random.randint(0,sal.shape[1],n) y = np.random.randint(0,sal.shape[0],n) z = np.random.rand(n)*max
これは結果の確認のためのコード無くても問題ない．画像uniは後で使わない．点かなり多め．．．
uni = np.zeros(shape = sal.shape, dtype="int8")
uni[y,x] = 1
plt.imshow(uni, cmap="gray")
plt.axis('off');
plt.subplots_adjust(left=0, right=1, bottom=0, top=1)
cv2.imwrite("uniform.png", rnd)
確率分布関数以上の点を削除
この記事の一番肝のコードがこれ． 1行目でrndには画像と同じサイズの中身が0の画像が作製される．int8型にしたのは後でor演算して画像に収めるため． 2行目では，$z < f(x,y) $を満たすとき1，満たさないとき0を，x,yのデータ数（ベクトルの次元数）分だけpntに代入している． なのでpntはxやyと同じshapeになる． 3行目では，画像rndの各画素$(x,y)$にpntの1か0が代入される．
rnd = np.zeros(shape = sal.shape, dtype="int8") pnt = np.where(z < sal[y,x], 1, 0) rnd[y,x] = rnd[y,x] | pnt
左が結果．なんか点が少ない．確かに$z_i < f(x_i,y_i) $を満たす点だけに絞られるので数がかなり減る．$z_{\rm max}$付近の値をとるのは$(x,y)$空間中のごく一部だけ．その他の点は殆ど無駄になる．これがこのアルゴリズムの難点． nを100倍のn = 10000000にして再実行すると次が得られる．それっぽい結果（右）が得られた．もっと効率の良い方法があると思う．
