Digital Reactor
統計・確率

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

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

はじめに

緯度と経度をそれぞれ一様乱数で決めて球面に点を置くと、点は両極に密集し、赤道付近はまばらになります。分子シミュレーションの初期配置や3Dグラフィックスのパーティクル配置でこの偏りを放置すると、極の近くだけ計算や描画が過密になります。均一に配置するにはひと工夫が要り、やり方は大きく2つあります。緯線の長さに合わせた乱数で緯度を決める方法と、黄金比を使って規則的に並べる方法です。どちらも短いPythonコードで書けます。

対象読者:

  • プログラミングの基礎知識（Python）がある方
  • シミュレーションや3Dグラフィックスに興味がある方
  • 数学的な概念（確率密度関数、累積分布関数、逆関数など）に抵抗がない方

記事のポイント:

  • 緯線の長さに比例した確率密度で緯度を生成し、点を均一に配置する方法を解説
  • 黄金比（フィボナッチ数列）を利用した、より簡便な点の配置方法を紹介
  • それぞれの方法をPythonコードで実装し、結果を可視化

緯線の長さに比例した確率で緯度を選ぶ

緯度と経度を単純な一様乱数で決めると点が偏るのは、緯線の長さが緯度によって違うためです。緯線は赤道で最も長く、極に近づくほど短くなるのに、どの緯度も同じ頻度で選ばれるので、短い緯線の上に同じ数の点が押し込まれます。均一な分布にするには、緯線の長さに比例した頻度で緯度を生成します。

緯線の長さから累積分布関数を作る

話を簡単にするため、半径1の単位球面を考えます。緯度をtt、経度をuuとします。ttは[−π/2,π/2][-\pi/2, \pi/2]、uuは[−π,π][-\pi, \pi]の範囲の値をとります。このとき、緯度ttにおける緯線の長さはcos⁡(t)\cos(t)で表されます。

緯線の長さに比例した確率で点を生成するため、累積分布関数f(x)f(x)を考えます。緯度−π/2-\pi/2からxxまでcos⁡(t)\cos(t)を積分し、全体で1になるよう規格化すると、次のようになります。

f(x)=12sin⁡(x)+12f(x) = \frac{1}{2} \sin(x) + \frac{1}{2}

2つの12\frac{1}{2}は規格化のための定数で、これによりf(−π/2)=0f(-\pi/2) = 0、f(π/2)=1f(\pi/2) = 1となり、ttの範囲内で確率の総和が1になります。

この累積分布関数をグラフで表すと、次のようになります。

累積分布関数の逆関数

次に、この累積分布関数の逆関数f−1(x)f^{-1}(x)を求めます。これにより、一様乱数から緯度ttを生成できます。逆関数は次のように表されます。

f−1(x)={−sin⁡−1(1−2x)(0<x<1)未定義(x≤0,x≥1)f^{-1}(x) = \begin{cases} - \sin^{-1}(1 - 2x) & (0 < x < 1) \\ \text{未定義} & (x \leq 0, x \geq 1) \end{cases}

ここで、sin⁡−1\sin^{-1}は逆正弦関数（アークサイン）を表します。xxが0から1の範囲の一様乱数であることに注意してください。

この逆関数をグラフで表すと、次のようになります。

乱数を用いた点群の生成と可視化 (Pythonコード)

上記の逆関数を用いて緯度ttを生成し、経度uuは[−π,π][-\pi, \pi]の一様乱数とすることで、球面上の点を生成します。以下にPythonコードを示します。

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

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()

このコードを実行すると、球面上に点が均一に分布している様子が確認できます。

黄金比を用いた球面上への点の配置

乱数を使わずに、黄金比を利用して球面上に点を配置する方法もあります。この方法は、ヒマワリの種の配列など、自然界に見られるフィボナッチ数列に基づいた配置方法です。

黄金比とフィボナッチ数列

黄金比は、(1+5)/2(1 + \sqrt{5}) / 2 で表される無理数で、約1.618です。フィボナッチ数列は、隣り合う2つの数の和が次の数になる数列（1, 1, 2, 3, 5, 8…）です。黄金比とフィボナッチ数列は密接な関係があり、この性質を利用して点を配置します。

Pythonコードによる実装

以下に、黄金比を用いて球面状に点を配置するPythonコードを示します。

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()

このコードを実行すると、点が球面全体に比較的均一に配置されていることがわかります。

まとめ

乱数を用いた方法と黄金比を用いた方法の違いを整理すると、次のようになります。

手法の概要特徴
乱数を用いた方法緯線の長さを考慮した確率密度関数を使用。逆関数を求める必要がある。点同士の距離はランダム。点が少ないと偏ることがある。
黄金比を用いた方法フィボナッチ数列に基づく配置。比較的シンプルな計算で実現できる。点同士の距離は均一。点が少なくても偏らない。

これらの方法は、分子動力学シミュレーションでの原子の初期配置、3Dモデリングでのパーティクルやオブジェクトの配置、機械学習でのデータのサンプリング、コンピュータグラフィックスでのテクスチャマッピングなど、球面上に点を散らす場面で広く使えます。点数が少ない場合や均一さを重視する場合は黄金比を用いた方法が扱いやすく、点同士がランダムに散らばってほしい場合は乱数を用いた方法が向きます。

関連記事

← 技術ブログ一覧へ