方針
簡単のために単位球を考えます。乱数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()
均一に点を播くことができました。何かシミュレーションをやるときの初期値として使えそうです。
黄金比を使った点の撒き方
コード
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()
