球面上に均一にランダム分布する点群を生成する方法を考えてみます。
方針
簡単のために単位球を考えます。乱数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} 1-\sqrt{1-x} & (0<x<1) \\ \text{未定義} & (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()
均一に点を播くことができました。何かシミュレーションをやるときの初期値として使えそうです。
黄金比を使った点の撒き方
こちらの方法では、ランダムではなく、点と点の間隔が均等に近い配置になります。
コード
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
def fib(N):
f = (np.sqrt(5)-1)/2
arr = np.linspace(-N, N, N*2+1)
theta = np.arcsin(arr/N)
phi = 2*np.pi*arr*f
x = np.cos(theta)*np.cos(phi)
y = np.cos(theta)*np.sin(phi)
z = np.sin(theta)
return x, y, z
x, y, z = fib(10000)
fig = plt.figure(figsize=(15,15))
ax = Axes3D(fig)
ax.plot(x, y, z, "o", ms=2, mew=0.5)
plt.show()