获取球体上均匀分布点的区域
Obtaining a region of evenly distributed points on a sphere
这个网站上有几个关于在球体表面分布点的问题，但所有这些都是基于实际生成该球体上的所有点。到目前为止，我最喜欢的是 Evenly distributing n points on a sphere.
中讨论的黄金螺旋
我需要用数万亿个点覆盖一个球体，但只需要实际生成一个很小的表面区域（地球低至约 10 米，观察大约 1 km^2 的区域）。为该区域生成的点必须与为整个球体生成的点匹配（即，将小区域拼接在一起必须产生与生成更大区域相同的结果），并且生成应该非常快。
我尝试使用具有如此大量点的黄金螺旋线的尝试因浮点精度问题而受阻。
我想出的最好办法是在等距纬度上生成点，并根据该纬度的圆周计算经度间距。然而，结果远非令人满意（尤其是由此产生的水平点环）。
有没有人有关于在球体表面生成一小块分布点的建议？
geodesic sphere 的顶点在此应用程序中可以很好地工作。
您从一个二十面体开始，将每个面分成任意分辨率的三角形网格，然后将这些点投影到球体表面。
Fibonacci sphere approximation 很容易有效地推广到点计算的子集，因为解析公式非常 straight-forward.
下面的代码在我较弱的笔记本电脑和相对优化不足的 python 实现上，在几秒钟的运行时间内计算了下面显示的一万亿个点的子集。
计算上面的代码如下，包括一种验证子集计算与brute-force计算完全相同的方法（但是不要尝试万亿点，它永远不会完成，除非你有一个 super-computer!)
请注意，当您对超过 10 亿个点进行计算时，使用 128 位双精度数是绝对要求，否则会出现主要的量化伪像！
运行时按 r' * N 缩放，其中 r' 是子集与整个球体的比率。因此，可以非常有效地计算出非常小的 r'。
#!/usr/bin/env python3
import argparse
import mpl_toolkits.mplot3d.axes3d as ax3d
import matplotlib.pyplot as plt
import numpy as np
def fibonacci_sphere_pts(num_pts):
ga = (3 - np.sqrt(5)) * np.pi # golden angle
# Create a list of golden angle increments along tha range of number of points
theta = ga * np.arange(num_pts)
# Z is a split into a range of -1 to 1 in order to create a unit circle
z = np.linspace(1 / num_pts - 1, 1 - 1 / num_pts, num_pts)
# a list of the radii at each height step of the unit circle
radius = np.sqrt(1 - z * z)
# Determine where xy fall on the sphere, given the azimuthal and polar angles
y = radius * np.sin(theta)
x = radius * np.cos(theta)
return np.asarray(list(zip(x,y,z)))
def fibonacci_sphere(num_pts):
x,y,z = zip(*fibonacci_sphere_subset(num_pts))
# Display points in a scatter plot
fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")
ax.scatter(x, y, z)
plt.show()
def fibonacci_sphere_subset_pts(num_pts, p0, r0 ):
"""
Get a subset of a full fibonacci_sphere
"""
ga = (3 - np.sqrt(5)) * np.pi # golden angle
x0, y0, z0 = p0
z_s = 1 / num_pts - 1
z_e = 1 - 1 / num_pts
# linspace formula for range [z_s,z_e] for N points is
# z_k = z_s + (z_e - z_s) / (N-1) * k , for k [0,N)
# therefore k = (z_k - z_s)*(N-1) / (z_e - z_s)
# would be the closest value of k
k = int(np.round((z0 - z_s) * (num_pts - 1) / (z_e - z_s)))
# here a sufficient number of "layers" of the fibonacci sphere must be
# selected to obtain enough points to be a superset of the subset given the
# radius, we use a heuristic to determine the number but it can be obtained
# exactly by the correct formula instead (by choosing an upperbound)
dz = (z_e - z_s) / (num_pts-1)
n_dk = int(np.ceil( r0 / dz ))
dk = np.arange(k - n_dk, k + n_dk+1)
dk = dk[np.where((dk>=0)&(dk<num_pts))[0]]
# NOTE: *must* use long double over regular doubles below, otherwise there
# are major quantization errors in the output for large number of points
theta = ga * dk.astype(np.longdouble)
z = z_s + (z_e - z_s ) / (num_pts-1) *dk
radius = np.sqrt(1 - z * z)
y = radius * np.sin(theta)
x = radius * np.cos(theta)
idx = np.where(np.square(x - x0) + np.square(y-y0) + np.square(z-z0) <= r0*r0)[0]
return x[idx],y[idx],z[idx]
def fibonacci_sphere_subset(num_pts, p0, r0, do_compare=False ):
"""
Display fib sphere subset points and optionally compare against bruteforce computation
"""
x,y,z = fibonacci_sphere_subset_pts(num_pts,p0,r0)
if do_compare:
subset = zip(x,y,z)
subset_bf = fibonacci_sphere_pts(num_pts)
x0,y0,z0 = p0
subset_bf = [ (x,y,z) for (x,y,z) in subset_bf if np.square(x - x0) + np.square(y-y0) + np.square(z-z0) <= r0*r0 ]
subset_bf = np.asarray(subset_bf)
if np.allclose(subset,subset_bf):
print('PASS: subset and bruteforce computation agree completely')
else:
print('FAIL: subset and bruteforce computation DO NOT agree completely')
# Display points in a scatter plot
fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")
ax.scatter(x, y, z)
plt.show()
if __name__ == "__main__":
parser = argparse.ArgumentParser(description="fibonacci sphere")
parser.add_argument(
"numpts", type=int, help="number of points to distribute along sphere"
)
args = parser.parse_args()
# hard-coded point to query with a tiny fixed radius
p0 = (.5,.5,np.sqrt(1. - .5*.5 - .5*.5)) # coordinate of query point representing center of subset, note all coordinates fall between -1 and 1
r0 = .00001 # the radius of the subset, a very small number is chosen as radius of full sphere is 1.0
fibonacci_sphere_subset(int(args.numpts),p0,r0,do_compare=False)
这个网站上有几个关于在球体表面分布点的问题，但所有这些都是基于实际生成该球体上的所有点。到目前为止，我最喜欢的是 Evenly distributing n points on a sphere.
中讨论的黄金螺旋我需要用数万亿个点覆盖一个球体，但只需要实际生成一个很小的表面区域（地球低至约 10 米，观察大约 1 km^2 的区域）。为该区域生成的点必须与为整个球体生成的点匹配（即，将小区域拼接在一起必须产生与生成更大区域相同的结果），并且生成应该非常快。
我尝试使用具有如此大量点的黄金螺旋线的尝试因浮点精度问题而受阻。
我想出的最好办法是在等距纬度上生成点，并根据该纬度的圆周计算经度间距。然而，结果远非令人满意（尤其是由此产生的水平点环）。
有没有人有关于在球体表面生成一小块分布点的建议？
geodesic sphere 的顶点在此应用程序中可以很好地工作。
您从一个二十面体开始，将每个面分成任意分辨率的三角形网格，然后将这些点投影到球体表面。
Fibonacci sphere approximation 很容易有效地推广到点计算的子集，因为解析公式非常 straight-forward.
下面的代码在我较弱的笔记本电脑和相对优化不足的 python 实现上，在几秒钟的运行时间内计算了下面显示的一万亿个点的子集。
计算上面的代码如下，包括一种验证子集计算与brute-force计算完全相同的方法（但是不要尝试万亿点，它永远不会完成，除非你有一个 super-computer!)
请注意，当您对超过 10 亿个点进行计算时，使用 128 位双精度数是绝对要求，否则会出现主要的量化伪像！
运行时按 r' * N 缩放，其中 r' 是子集与整个球体的比率。因此，可以非常有效地计算出非常小的 r'。
#!/usr/bin/env python3
import argparse
import mpl_toolkits.mplot3d.axes3d as ax3d
import matplotlib.pyplot as plt
import numpy as np
def fibonacci_sphere_pts(num_pts):
ga = (3 - np.sqrt(5)) * np.pi # golden angle
# Create a list of golden angle increments along tha range of number of points
theta = ga * np.arange(num_pts)
# Z is a split into a range of -1 to 1 in order to create a unit circle
z = np.linspace(1 / num_pts - 1, 1 - 1 / num_pts, num_pts)
# a list of the radii at each height step of the unit circle
radius = np.sqrt(1 - z * z)
# Determine where xy fall on the sphere, given the azimuthal and polar angles
y = radius * np.sin(theta)
x = radius * np.cos(theta)
return np.asarray(list(zip(x,y,z)))
def fibonacci_sphere(num_pts):
x,y,z = zip(*fibonacci_sphere_subset(num_pts))
# Display points in a scatter plot
fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")
ax.scatter(x, y, z)
plt.show()
def fibonacci_sphere_subset_pts(num_pts, p0, r0 ):
"""
Get a subset of a full fibonacci_sphere
"""
ga = (3 - np.sqrt(5)) * np.pi # golden angle
x0, y0, z0 = p0
z_s = 1 / num_pts - 1
z_e = 1 - 1 / num_pts
# linspace formula for range [z_s,z_e] for N points is
# z_k = z_s + (z_e - z_s) / (N-1) * k , for k [0,N)
# therefore k = (z_k - z_s)*(N-1) / (z_e - z_s)
# would be the closest value of k
k = int(np.round((z0 - z_s) * (num_pts - 1) / (z_e - z_s)))
# here a sufficient number of "layers" of the fibonacci sphere must be
# selected to obtain enough points to be a superset of the subset given the
# radius, we use a heuristic to determine the number but it can be obtained
# exactly by the correct formula instead (by choosing an upperbound)
dz = (z_e - z_s) / (num_pts-1)
n_dk = int(np.ceil( r0 / dz ))
dk = np.arange(k - n_dk, k + n_dk+1)
dk = dk[np.where((dk>=0)&(dk<num_pts))[0]]
# NOTE: *must* use long double over regular doubles below, otherwise there
# are major quantization errors in the output for large number of points
theta = ga * dk.astype(np.longdouble)
z = z_s + (z_e - z_s ) / (num_pts-1) *dk
radius = np.sqrt(1 - z * z)
y = radius * np.sin(theta)
x = radius * np.cos(theta)
idx = np.where(np.square(x - x0) + np.square(y-y0) + np.square(z-z0) <= r0*r0)[0]
return x[idx],y[idx],z[idx]
def fibonacci_sphere_subset(num_pts, p0, r0, do_compare=False ):
"""
Display fib sphere subset points and optionally compare against bruteforce computation
"""
x,y,z = fibonacci_sphere_subset_pts(num_pts,p0,r0)
if do_compare:
subset = zip(x,y,z)
subset_bf = fibonacci_sphere_pts(num_pts)
x0,y0,z0 = p0
subset_bf = [ (x,y,z) for (x,y,z) in subset_bf if np.square(x - x0) + np.square(y-y0) + np.square(z-z0) <= r0*r0 ]
subset_bf = np.asarray(subset_bf)
if np.allclose(subset,subset_bf):
print('PASS: subset and bruteforce computation agree completely')
else:
print('FAIL: subset and bruteforce computation DO NOT agree completely')
# Display points in a scatter plot
fig = plt.figure()
ax = fig.add_subplot(111, projection="3d")
ax.scatter(x, y, z)
plt.show()
if __name__ == "__main__":
parser = argparse.ArgumentParser(description="fibonacci sphere")
parser.add_argument(
"numpts", type=int, help="number of points to distribute along sphere"
)
args = parser.parse_args()
# hard-coded point to query with a tiny fixed radius
p0 = (.5,.5,np.sqrt(1. - .5*.5 - .5*.5)) # coordinate of query point representing center of subset, note all coordinates fall between -1 and 1
r0 = .00001 # the radius of the subset, a very small number is chosen as radius of full sphere is 1.0
fibonacci_sphere_subset(int(args.numpts),p0,r0,do_compare=False)