Grids on a sphere

球面上で定義される関数を離散化するためのグリッドを定義する方法を考えます. ここでは \(N\) 個の点を球面上に配置することを考えます. ここで, 隣接する点の距離ができるだけ均一になるように配置することを考えます. こうした点列を生成するためのアルゴリズムが Stack Overflow にて議論されていました. このページでは議論されていた手法について日本で解説します.

Fibonacci Sphere

このアルゴリズムは González (2010) にて詳しく解説されています. このアルゴリズムでは球面上に 2 種類の螺旋を描き, その交点を点列として採用することで球面をカバーします. まずはデモンストレーションとして 2 次元平面での単位円を埋め尽くす点列を生成してみます. Python での実装例を以下に示します.

import numpy as np

n_point = 100
index = np.arange(n_point) + 0.5
r = np.sqrt(index / n_point)

φ = (np.sqrt(5) + 1) / 2.0
θ = 2 * np.pi / φ * index

x = r * np.cos(θ)
y = r * np.sin(θ)

以下に \(N = 20, 1000\) で生成した点列を示します.

点列は 2 本の螺旋の交点で定義されています. \(\phi\) を黄金比とすると 1 本目の螺旋 (青色) は以下の式で定義されています.

\[ r = \sqrt{s \over N}, \quad \theta = 2\pi \phi^{-1} s \]

また, 2 本目の螺旋 (オレンジ色) は以下の式で定義されています.

\[ r = \sqrt{t \over N}, \quad \theta = \pi - 2\pi \phi^{-2} t \]

ここで \(s, t\) は媒介変数です. 2 本の螺旋の交点を求めると以下の式を得ます.

\[ s = t = 0.5 + n \quad (n = 0, 1, \ldots) \]

\(s = N - 0.5\) で止めれば半径 1 の円盤内にとどまりますが, 螺旋はいくらでも広げることができます. ここで \(r \approx \sqrt{s}\) かつ \(\theta \propto s\) なので, 半径 \(r\) の円盤内に含まれる交点の数は円盤の面積に比例します. つまり, このように生成した点列は 2 次元平面内をほぼ一様に埋め尽くすことができます.

ここで \(N\) を大きくすると点列全体が原点に向かって相似収縮します. 上図で見たように単位円盤内を効率よく埋め尽くす点列を生成することができます. このパターンはヒマワリの種の配置に似ていることから sunflower spiral とも呼ばれます.

このパターンを 3 次元の球面に展開することを考えます. 上半面と下半面で 2 つの sunflower spirals を用意して貼り合わせることで対応できそうです. 球面上で点列を生成するアルゴリズムを以下に示します.

import numpy as np

def fibonacci_sphere(n_side):
    point = []
    n_point = 2 * n_side + 1
    φ = (np.sqrt(5) + 1) / 2.0

    for i in range(-n_side, n_side + 1):
        z = i / n_point
        radius = np.sqrt(1 - z * z)  # radius at z

        θ = 2 * np.pi / φ * i  # golden angle increment

        x = np.cos(θ) * radius
        y = np.sin(θ) * radius

        point.append((x, y, z))

    return np.array(point)

\(N = 2n + 1, i = {-n, -n+1, \ldots, n}\) として, 点列を \(z\) 座標が大きい順に並び替えると, 各点の緯度と経度は以下のように計算できます. \(z\) 座標については等間隔に並ぶこと, また \(z\) 軸の極には点が配置されないと言った特徴があります.

\[ \mathrm{lat}_i = \arcsin\left(1 - \frac{i}{2n + 1} \right), \quad \mathrm{lon}_i = 2\pi i \phi^{-1}. \]

この点列を球面上にプロットしたものを以下に示します.

各点の最近傍点までの距離の頻度分布を以下に示します. n_side が大きくなると最近傍点までの距離はほぼ反比例して減少します. 距離の標準偏差を平均値で規格化した値が \(\Delta\) です. この値が小さいほど点列の分布が均一であることを示します. n_side が大きくなると \(\Delta\) が若干大きくなる傾向はありますが, おおむね 1 % 以下に抑えられていることがわかります. Fibonacci Sphere は極めて単純なアルゴリズムで球面を埋め尽くす点列を生成することができます.

この Grid を使用して Gaussian Process で閉じたポリゴンを生成してみます.

import numpy as np

p = fibonacci_sphere(100)

mu = 2.0 * np.ones(p.shape[0])
cov = np.exp(0.8 * np.dot(p, p.T)) + np.eye(p.shape[0]) * 1e-2

q = p * np.random.multivariate_normal(mu, cov)[:, None]

p が生成した grid 点, q が Gaussian Process によって生成された点列です. カーネルには dot product (von Mises-Fisher) kernel を採用しました. 以下では生成した点列を 3D ポリゴンとして可視化しています.