球面上に均一にランダム分布する点群を生成する

方針

簡単のために単位球を考えます。乱数tで緯度を決め、乱数uで経度を決めます。ただし、緯度ごとに、緯線の長さが異なるため、球面上に均一に分布させるために、この緯線の長さに比例した頻度でtを生成する必要があります。

なお、乱数を使わない方法としては、フィボナッチ数列を使った配置方法があります(ひまわりの種の配置と同じです)。

緯線の長さに比例した確率密度関数をもつ乱数の生成

tの範囲は[-pi/2, pi/2]、uの範囲は[-pi, pi]とします。緯線の長さは、cos(t)で計算されます。

t=pi/2を起点とする累積分布関数f(x)は、これの積分により、f(x) = sin(x)*1/2 + 1/2 となります。2つの1/2は、tの範囲内で、確率分布の合計を1とし、f(-pi/2)=0, f(pi/2)=1を満たすようにするための係数です。

プロットすると下図のようになります。

ここで、このページの内容を使います。f(x)の逆関数を求めます。

一般解, (nを整数として)
$$ f^{-1}(x) = 2\pi n – sin^{-1}(1 – 2x) $$

ですが、0 < x < 1 に注意すると、
$$ f^{-1}(x) =
\begin{cases}
– sin^{-1}(1 – 2x) &(0<x<1) \\
未定義 &(x \leq 0, x \geq 1)
\end{cases}
$$
となります。

点群の生成と可視化

最後に確認のために、このようにして得られた変換式をtに適用し、uは一様乱数のままで、球面上に点を生成し、可視化してみます。このページの最初に表示している図が出力されます。

コード

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import random

points = []
for i in range(10000):
    t = random.random()
    t = - np.arcsin(1-2*t)
    u = random.random() * 2 * np.pi - np.pi

    x = np.cos(t) * np.cos(u)
    y = np.cos(t) * np.sin(u)
    z = np.sin(t)
    points.append([x, y, z])
points = np.array(points)

fig = plt.figure(figsize=(15,15))
ax = Axes3D(fig)
ax.plot(points[:, 0], points[:,1], points[:, 2], "o", ms=2, mew=0.5)
plt.show()

均一に点を播くことができました。何かシミュレーションをやるときの初期値として使えそうです。