如何使用 Python 从 3-D 椭球生成点的随机样本？
How to generate a random sample of points from a 3-D ellipsoid using Python?
我正在尝试从 3-D 椭球体中均匀地采样大约 1000 个点。有什么方法可以编码，这样我们就可以从椭球方程开始得到点吗？
我想要椭圆体表面的点。
考虑使用 Monte-Carlo 模拟：生成随机 3D 点；检查点是否在椭圆体内；如果是，请保留它。重复直到获得 1,000 分。
P.S。由于 OP 更改了他们的问题，因此此答案不再有效。
J.F。 Williamson，“Random selection of points distributed on curve surfaces”，Physics in Medicine & Biology 32(10)，1987，描述了在参数化表面上选择均匀随机点的一般方法。它是一种 acceptance/rejection 方法，根据其拉伸因子 (norm-of-gradient) 接受或拒绝每个候选点。要将此方法用于参数化曲面，必须了解有关曲面的几件事，即 -
x(u, v)、y(u, v)和z(u, v)，它们是从二维坐标u和[=14生成3维坐标的函数=],
范围u和v,
g(point)，曲面上每个点的梯度范数（“拉伸因子”），以及
gmax，整个面的最大值g。
那么算法就是：
- 在曲面上生成一个点，
xyz。
- 如果
g(xyz) >= RNDU01()*gmax，其中RNDU01()是[0, 1)中的均匀随机变量，接受该点。否则，重复此过程。
Chen 和 Glotzer (2007) 在“延长病毒衣壳形成的现象学模型的模拟研究”中将该方法应用于长球体（椭圆体的一种形式）的表面，物理评论 E 75(5), 051504 (preprint).
这是一个通用函数，用于在球体、椭球体或任何具有 a、b 和 c 参数的三轴椭球体的表面上选取一个随机点。请注意，直接生成角度不会提供均匀分布，并且会导致沿 z 方向的点过多。相反，phi 是作为随机生成的 cos(phi) 的倒数获得的。
import numpy as np
def random_point_ellipsoid(a,b,c):
u = np.random.rand()
v = np.random.rand()
theta = u * 2.0 * np.pi
phi = np.arccos(2.0 * v - 1.0)
sinTheta = np.sin(theta);
cosTheta = np.cos(theta);
sinPhi = np.sin(phi);
cosPhi = np.cos(phi);
rx = a * sinPhi * cosTheta;
ry = b * sinPhi * sinTheta;
rz = c * cosPhi;
return rx, ry, rz
这个函数是从这个post:https://karthikkaranth.me/blog/generating-random-points-in-a-sphere/
理论
使用this excellent answer to the MSE question How to generate points uniformly distributed on the surface of an ellipsoid?我们可以
generate a point uniformly on the sphere, apply the mapping f :
(x,y,z) -> (x'=ax,y'=by,z'=cz) and then correct the distortion
created by the map by discarding the point randomly with some
probability p(x,y,z).
假设椭圆体的 3 个轴被命名为
0 < a < b < c
我们用
丢弃生成的点
p(x,y,z) = 1 - mu(x,y,y)/mu_max
概率，即我们用mu(x,y,y)/mu_max概率保持它
mu(x,y,z) = ((acy)^2 + (abz)^2 + (bcx)^2)^0.5
和
mu_max = bc
实施
import numpy as np
np.random.seed(42)
# Function to generate a random point on a uniform sphere
# (relying on
def randompoint(ndim=3):
vec = np.random.randn(ndim,1)
vec /= np.linalg.norm(vec, axis=0)
return vec
# Give the length of each axis (example values):
a, b, c = 1, 2, 4
# Function to scale up generated points using the function `f` mentioned above:
f = lambda x,y,z : np.multiply(np.array([a,b,c]),np.array([x,y,z]))
# Keep the point with probability `mu(x,y,z)/mu_max`, ie
def keep(x, y, z, a=a, b=b, c=c):
mu_xyz = ((a * c * y) ** 2 + (a * b * z) ** 2 + (b * c * x) ** 2) ** 0.5
return mu_xyz / (b * c) > np.random.uniform(low=0.0, high=1.0)
# Generate points until we have, let's say, 1000 points:
n = 1000
points = []
while len(points) < n:
[x], [y], [z] = randompoint()
if keep(x, y, z):
points.append(f(x, y, z))
支票
检查生成的所有点是否满足椭球条件（即x^2/a^2 + y^2/b^2 + z^2/c^2 = 1）：
for p in points:
pscaled = np.multiply(p,np.array([1/a,1/b,1/c]))
assert np.allclose(np.sum(np.dot(pscaled,pscaled)),1)
运行时没有出现任何错误。可视化点：
import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.add_subplot(projection="3d")
points = np.array(points)
ax.scatter(points[:, 0], points[:, 1], points[:, 2])
# set aspect ratio for the axes using
ax.set_box_aspect((np.ptp(points[:, 0]), np.ptp(points[:, 1]), np.ptp(points[:, 2])))
plt.show()
这些点似乎分布均匀。
有问题
在球体上生成一个点，然后将其重新投影而不对椭圆进行任何进一步校正，将导致分布失真。这与将 this posts's p(x,y,z) 设置为 0 基本相同。想象一个椭圆体，其中一个轴比另一个轴大几个数量级。这样，很容易看出，朴素的重投影是行不通的。
为了在椭球表面均匀生成点（“均匀”指的是面积密度），我们可以考虑使用球坐标的参数方程（来自Wikipedia）：
其中s = 1指的是semi-axesa, b, c给出的椭球体。
面积密度由体积元素给出
因此与 sin(theta) 成正比（对于 0 <= theta <= pi）。我们可以以此作为概率分布的基础，表示对于给定的theta值应该生成“多少”个点：其中面积密度为low/high，生成相应值的概率theta 也应该是 low/high。
因此，我们可以使用函数f(theta) = sin(theta)/2作为我们在区间[0, pi]上的概率分布。对应的累积分布函数为F(theta) = (1 - cos(theta))/2。现在我们可以使用 Inverse transform sampling 根据 f(theta) 从均匀随机分布生成 theta 的值。 phi的值可以直接从[0, 2*pi].
上的均匀分布得到
示例代码：
import matplotlib.pyplot as plt
import numpy as np
from numpy import sin, cos, pi
rng = np.random.default_rng(seed=0)
a, b, c = 10, 3, 1
N = 5000
phi = rng.uniform(0, 2*pi, size=N)
theta = np.arccos(1 - 2*rng.random(size=N))
x = a*sin(theta)*cos(phi)
y = b*sin(theta)*sin(phi)
z = c*cos(theta)
fig = plt.figure()
ax = fig.add_subplot(projection='3d')
ax.scatter(x, y, z, s=2)
plt.show()
产生以下情节：
一种适用于任何形状或表面的方法是将表面转换为任意高分辨率的体素表示（分辨率越高越好，但也越慢）。然后，您可以根据需要轻松 select 随机体素，然后您可以使用参数方程 select 体素内表面上的一个点。体素 selection 应该是完全无偏的，并且体素内点的 selection 将遭受与使用参数方程相同的偏差，但是如果有足够的体素，那么这些的大小偏差会很小。
您需要高质量的立方体相交代码，但需要像椭圆体这样的东西可以很容易地优化。我建议逐步通过细分为体素的边界框。快速距离检查将消除大多数立方体，您可以对可能存在交叉的立方体进行适当的交叉检查。对于立方体中的点，我很想做一些简单的事情，比如距中心的随机 XYZ 距离，然后从椭球体的中心投射一条光线，selected 点是光线与表面相交的地方.正如我上面所说，它会有偏差，但对于小体素，偏差可能足够小。
有些库可以非常有效地进行凸形相交，cube/elipsoid 将是其中一种选择。它们将被高度优化，但我认为距离剔除可能值得手动完成。您将需要一个库来区分表面相交和一个对象完全位于另一个对象内部。
如果你知道你的椭圆体与轴对齐，那么你可以很容易地做 voxel/edge 相交作为一堆 2D 正方形相交椭圆问题，要测试的正方形集定义为那些与上层中的那些相邻。那可能会更快。
使这种方法更易于管理的一个原因是您不需要为边缘情况编写所有代码（要解决可能导致丢失或丢失的浮点不准确问题需要大量工作在交叉点处加倍体素）。那是因为这些会非常罕见，所以它们不会影响您的抽样。
简单地找到椭圆内的所有体素然后丢弃具有 6 个邻居的所有体素甚至可能更快......有很多选择。这完全取决于性能的重要性。这将比其他建议慢得多，但如果您想要 ~1000 个点，那么 ~100,000 体素感觉是表面的最小值，因此您的边界框可能需要 ~1,000,000 体素。然而，即使在现代计算机上测试 1,000,000 个交叉点也相当快。
我正在尝试从 3-D 椭球体中均匀地采样大约 1000 个点。有什么方法可以编码，这样我们就可以从椭球方程开始得到点吗？
我想要椭圆体表面的点。
考虑使用 Monte-Carlo 模拟：生成随机 3D 点；检查点是否在椭圆体内；如果是，请保留它。重复直到获得 1,000 分。
P.S。由于 OP 更改了他们的问题，因此此答案不再有效。
J.F。 Williamson，“Random selection of points distributed on curve surfaces”，Physics in Medicine & Biology 32(10)，1987，描述了在参数化表面上选择均匀随机点的一般方法。它是一种 acceptance/rejection 方法，根据其拉伸因子 (norm-of-gradient) 接受或拒绝每个候选点。要将此方法用于参数化曲面，必须了解有关曲面的几件事，即 -
x(u, v)、y(u, v)和z(u, v)，它们是从二维坐标u和[=14生成3维坐标的函数=],范围
u和v,g(point)，曲面上每个点的梯度范数（“拉伸因子”），以及gmax，整个面的最大值g。
那么算法就是：
- 在曲面上生成一个点，
xyz。 - 如果
g(xyz) >= RNDU01()*gmax，其中RNDU01()是[0, 1)中的均匀随机变量，接受该点。否则，重复此过程。
Chen 和 Glotzer (2007) 在“延长病毒衣壳形成的现象学模型的模拟研究”中将该方法应用于长球体（椭圆体的一种形式）的表面，物理评论 E 75(5), 051504 (preprint).
这是一个通用函数，用于在球体、椭球体或任何具有 a、b 和 c 参数的三轴椭球体的表面上选取一个随机点。请注意，直接生成角度不会提供均匀分布，并且会导致沿 z 方向的点过多。相反，phi 是作为随机生成的 cos(phi) 的倒数获得的。
import numpy as np
def random_point_ellipsoid(a,b,c):
u = np.random.rand()
v = np.random.rand()
theta = u * 2.0 * np.pi
phi = np.arccos(2.0 * v - 1.0)
sinTheta = np.sin(theta);
cosTheta = np.cos(theta);
sinPhi = np.sin(phi);
cosPhi = np.cos(phi);
rx = a * sinPhi * cosTheta;
ry = b * sinPhi * sinTheta;
rz = c * cosPhi;
return rx, ry, rz
这个函数是从这个post:https://karthikkaranth.me/blog/generating-random-points-in-a-sphere/
理论
使用this excellent answer to the MSE question How to generate points uniformly distributed on the surface of an ellipsoid?我们可以
generate a point uniformly on the sphere, apply the mapping f : (x,y,z) -> (x'=ax,y'=by,z'=cz) and then correct the distortion created by the map by discarding the point randomly with some probability p(x,y,z).
假设椭圆体的 3 个轴被命名为
0 < a < b < c
我们用
丢弃生成的点p(x,y,z) = 1 - mu(x,y,y)/mu_max
概率，即我们用mu(x,y,y)/mu_max概率保持它
mu(x,y,z) = ((acy)^2 + (abz)^2 + (bcx)^2)^0.5
和
mu_max = bc
实施
import numpy as np
np.random.seed(42)
# Function to generate a random point on a uniform sphere
# (relying on
def randompoint(ndim=3):
vec = np.random.randn(ndim,1)
vec /= np.linalg.norm(vec, axis=0)
return vec
# Give the length of each axis (example values):
a, b, c = 1, 2, 4
# Function to scale up generated points using the function `f` mentioned above:
f = lambda x,y,z : np.multiply(np.array([a,b,c]),np.array([x,y,z]))
# Keep the point with probability `mu(x,y,z)/mu_max`, ie
def keep(x, y, z, a=a, b=b, c=c):
mu_xyz = ((a * c * y) ** 2 + (a * b * z) ** 2 + (b * c * x) ** 2) ** 0.5
return mu_xyz / (b * c) > np.random.uniform(low=0.0, high=1.0)
# Generate points until we have, let's say, 1000 points:
n = 1000
points = []
while len(points) < n:
[x], [y], [z] = randompoint()
if keep(x, y, z):
points.append(f(x, y, z))
支票
检查生成的所有点是否满足椭球条件（即x^2/a^2 + y^2/b^2 + z^2/c^2 = 1）：
for p in points:
pscaled = np.multiply(p,np.array([1/a,1/b,1/c]))
assert np.allclose(np.sum(np.dot(pscaled,pscaled)),1)
运行时没有出现任何错误。可视化点：
import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.add_subplot(projection="3d")
points = np.array(points)
ax.scatter(points[:, 0], points[:, 1], points[:, 2])
# set aspect ratio for the axes using
ax.set_box_aspect((np.ptp(points[:, 0]), np.ptp(points[:, 1]), np.ptp(points[:, 2])))
plt.show()
这些点似乎分布均匀。
在球体上生成一个点，然后将其重新投影而不对椭圆进行任何进一步校正，将导致分布失真。这与将 this posts's p(x,y,z) 设置为 0 基本相同。想象一个椭圆体，其中一个轴比另一个轴大几个数量级。这样，很容易看出，朴素的重投影是行不通的。
为了在椭球表面均匀生成点（“均匀”指的是面积密度），我们可以考虑使用球坐标的参数方程（来自Wikipedia）：
其中s = 1指的是semi-axesa, b, c给出的椭球体。
面积密度由体积元素给出
因此与 sin(theta) 成正比（对于 0 <= theta <= pi）。我们可以以此作为概率分布的基础，表示对于给定的theta值应该生成“多少”个点：其中面积密度为low/high，生成相应值的概率theta 也应该是 low/high。
因此，我们可以使用函数f(theta) = sin(theta)/2作为我们在区间[0, pi]上的概率分布。对应的累积分布函数为F(theta) = (1 - cos(theta))/2。现在我们可以使用 Inverse transform sampling 根据 f(theta) 从均匀随机分布生成 theta 的值。 phi的值可以直接从[0, 2*pi].
示例代码：
import matplotlib.pyplot as plt
import numpy as np
from numpy import sin, cos, pi
rng = np.random.default_rng(seed=0)
a, b, c = 10, 3, 1
N = 5000
phi = rng.uniform(0, 2*pi, size=N)
theta = np.arccos(1 - 2*rng.random(size=N))
x = a*sin(theta)*cos(phi)
y = b*sin(theta)*sin(phi)
z = c*cos(theta)
fig = plt.figure()
ax = fig.add_subplot(projection='3d')
ax.scatter(x, y, z, s=2)
plt.show()
产生以下情节：
一种适用于任何形状或表面的方法是将表面转换为任意高分辨率的体素表示（分辨率越高越好，但也越慢）。然后，您可以根据需要轻松 select 随机体素，然后您可以使用参数方程 select 体素内表面上的一个点。体素 selection 应该是完全无偏的，并且体素内点的 selection 将遭受与使用参数方程相同的偏差，但是如果有足够的体素，那么这些的大小偏差会很小。
您需要高质量的立方体相交代码，但需要像椭圆体这样的东西可以很容易地优化。我建议逐步通过细分为体素的边界框。快速距离检查将消除大多数立方体，您可以对可能存在交叉的立方体进行适当的交叉检查。对于立方体中的点，我很想做一些简单的事情，比如距中心的随机 XYZ 距离，然后从椭球体的中心投射一条光线，selected 点是光线与表面相交的地方.正如我上面所说，它会有偏差，但对于小体素，偏差可能足够小。
有些库可以非常有效地进行凸形相交，cube/elipsoid 将是其中一种选择。它们将被高度优化，但我认为距离剔除可能值得手动完成。您将需要一个库来区分表面相交和一个对象完全位于另一个对象内部。
如果你知道你的椭圆体与轴对齐，那么你可以很容易地做 voxel/edge 相交作为一堆 2D 正方形相交椭圆问题，要测试的正方形集定义为那些与上层中的那些相邻。那可能会更快。
使这种方法更易于管理的一个原因是您不需要为边缘情况编写所有代码（要解决可能导致丢失或丢失的浮点不准确问题需要大量工作在交叉点处加倍体素）。那是因为这些会非常罕见，所以它们不会影响您的抽样。
简单地找到椭圆内的所有体素然后丢弃具有 6 个邻居的所有体素甚至可能更快......有很多选择。这完全取决于性能的重要性。这将比其他建议慢得多，但如果您想要 ~1000 个点，那么 ~100,000 体素感觉是表面的最小值，因此您的边界框可能需要 ~1,000,000 体素。然而，即使在现代计算机上测试 1,000,000 个交叉点也相当快。