我想从一个三维椭球体中均匀地采样约1000个点。是否有一种编程方法，可以通过椭球体的方程获得这些点？
我想要位于椭球体表面上的点。
6个回答
6
理论
使用这个很好的答案来回答如何在椭球表面上生成均匀分布的点？的MSE问题，我们可以：
在球体上均匀生成一个点，应用映射f: （x，y，z）-> (x'=ax，y'=by，z'=cz)，然后通过以一定的概率p(x,y,z)随机丢弃点来纠正地图创建的扭曲。
假设椭球的三个轴被命名为：
0 < a < b < c
我们会丢弃生成的点，如果：
p(x,y,z) = 1 - mu(x,y,y)/mu_max
概率，即我们以 mu(x,y,z)/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 https://dev59.com/G1sX5IYBdhLWcg3wUuCa#33977530)
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 https://dev59.com/pGsz5IYBdhLWcg3wDz0A#64453375
ax.set_box_aspect((np.ptp(points[:, 0]), np.ptp(points[:, 1]), np.ptp(points[:, 2])))
plt.show()
这些点看起来分布均匀。
目前被接受的答案存在问题
在球面上生成一个点，然后没有任何进一步的修正就将其重新投影到椭圆上会导致分布失真。这基本上与将此帖子中的p（x，y，z）设置为0是一样的。想象一个轴比另一个轴大几个数量级的椭球体。很容易看出，朴素的重投影方法不会奏效。
- zabop - we're hiring
1
2我相信这是J.F. Williamson在1987年的《医学与生物物理学中的物理》杂志上建议的方法，用于在曲面上分布的点的随机选择。 - Peter O.
3
考虑使用蒙特卡罗模拟：生成一个随机的三维点，检查该点是否在椭球内部，如果是，则保留该点。重复此过程直到得到1,000个点。
- DYZ
3
J.F. Williamson在1987年的文章“Random selection of points distributed on curved surfaces”（发表于《物理医学与生物学》32卷10期），描述了一种通用方法来选择参数曲面上均匀随机的点。这是一种接受/拒绝的方法，根据其拉伸因子（梯度范数）接受或拒绝每个候选点。要将此方法用于参数曲面，必须了解有关曲面的几件事，即：
x(u,v)、y(u,v)和z(u,v)，它们是从二维坐标u和v生成三维坐标的函数，u和v的范围，g(point)，曲面上每个点的梯度范数（“拉伸因子”），以及gmax，整个曲面的最大g值。
然后算法如下：
- 在表面上生成一个点，
xyz。 - 如果
g(xyz) >= RNDU01()*gmax，其中RNDU01()是 [0,1) 中的均匀随机变量，则接受该点。否则，请重复此过程。
陈和Glotzer（2007）在“Simulation studies of a phenomenological model for elongated virus capsid formation”，Physical Review E 75（5），051504（preprint）中将该方法应用于一个长椭球体（一种椭球体形式）的表面。
- Peter O.
1
这是一个通用函数，用于在球体、椭球体或任何三轴椭球体的表面上选择一个随机点，其中参数为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
这个函数是从这篇文章中采用的：https://karthikkaranth.me/blog/generating-random-points-in-a-sphere/
- Nikolay Frick
1
2这不就是在球体上创建随机点，然后重新投影到椭球体上吗？如果是的话，为什么最终分布会均匀并不明显。事实上，它并不均匀。 - zabop
0
根据“uniformly”所指的内容，不同的方法适用。无论如何，我们可以使用球坐标系下的参数方程（来自Wikipedia）：
其中s = 1表示由半轴a > b > c给定的椭球体。通过这些方程，我们可以推导出相关的体积/面积元素并生成点，使得它们被生成的概率与该体积/面积元素成比例。这将在椭球体表面提供恒定的体积/面积密度。
1. 恒定体积密度
此方法在椭球体表面生成点，使得它们在整个表面上的体积密度是恒定的。其结果是一维投影（即x、y、z坐标）均匀分布；详见下图。
三轴椭球体的体积元素为（参见此处）：
因此，它与sin(theta)成比例（对于0 <= theta <= pi）。我们可以将其作为概率分布的基础，指示针对给定值theta应生成“多少”点：在密度低/高的区域，生成相应值theta的概率也应该低/高。
因此，我们可以使用函数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()
它生成以下绘图：
x，y，z的密度图）：import seaborn as sns
sns.kdeplot(data=dict(x=x, y=y, z=z))
plt.show()
2. 常数面密度
该方法生成椭球体表面上的点，使它们的面积密度在整个椭球体表面上保持不变。
同样地，我们首先计算相应的面积元素。为了简单起见，我们可以使用SymPy：
from sympy import cos, sin, symbols, Matrix
a, b, c, t, p = symbols('a b c t p')
x = a*sin(t)*cos(p)
y = b*sin(t)*sin(p)
z = c*cos(t)
J = Matrix([
[x.diff(t), x.diff(p)],
[y.diff(t), y.diff(p)],
[z.diff(t), z.diff(p)],
])
print((J.T @ J).det().simplify())
这将产生
-a**2*b**2*sin(t)**4 + a**2*b**2*sin(t)**2 + a**2*c**2*sin(p)**2*sin(t)**4 - b**2*c**2*sin(p)**2*sin(t)**4 + b**2*c**2*sin(t)**4
并且进一步简化为（除以（a*b）**2并取sqrt）：
sin(t)*np.sqrt(1 + ((c/b)**2*sin(p)**2 + (c/a)**2*cos(p)**2 - 1)*sin(t)**2)
由于在这种情况下，区域元素比较复杂，我们可以使用拒绝抽样：
import matplotlib.pyplot as plt
import numpy as np
from numpy import cos, sin
def f_redo(t, p):
return (
sin(t)*np.sqrt(1 + ((c/b)**2*sin(p)**2 + (c/a)**2*cos(p)**2 - 1)*sin(t)**2)
< rng.random(size=t.size)
)
rng = np.random.default_rng(seed=0)
N = 5000
a, b, c = 10, 3, 1
t = rng.uniform(0, np.pi, size=N)
p = rng.uniform(0, 2*np.pi, size=N)
redo = f_redo(t, p)
while redo.any():
t[redo] = rng.uniform(0, np.pi, size=redo.sum())
p[redo] = rng.uniform(0, 2*np.pi, size=redo.sum())
redo[redo] = f_redo(t[redo], p[redo])
x = a*np.sin(t)*np.cos(p)
y = b*np.sin(t)*np.sin(p)
z = c*np.cos(t)
fig = plt.figure()
ax = fig.add_subplot(projection='3d')
ax.scatter(x, y, z, s=2)
plt.show()
该操作产生以下分布：
x，y，z）：
- a_guest
2
你能否为theta和phi的生成添加更多的解释？对我来说，theta和phi独立于a、b、c似乎是错误的。如果你考虑一个非常拉长的椭球体，很明显这不应该是这种情况。 - zabop
@zabop
theta和phi的生成基于体积元素dxdydz，该元素取决于a,b,c；然而，这只是一个比例关系，因此对于推导相应的概率函数并不相关。该方法基于所生成点的体积密度在椭球面上是恒定的；这是该方法的“均匀”特性，并且对于任何椭球都是有效的。也许您使用了不同的“均匀”概念？如果是这样，请更精确地说明您的方法中使分布均匀的确切因素是什么？ - a_guest0
一种通用的方法是将任何形状或表面转换为任意高分辨率的体素表示（分辨率越高越好，但速度也越慢）。然后您可以随意选择体素，然后使用参数方程在体素内选择表面上的点。体素选择应该完全无偏见，而在体素内选择点将遭受使用参数方程带来的同样偏差，但如果有足够的体素，则这些偏差的大小将非常小。
您需要一个高质量的立方体交集代码，但对于像椭球这样的形状，可以很容易地进行优化。我建议逐步遍历细分为体素的边界框。快速距离检查将消除大多数立方体，并且您可以对可能存在交集的立方体进行适当的交集检查。对于立方体内的点，我倾向于做一些简单的事情，例如从中心开始的随机XYZ距离，然后从椭球的中心投射一条射线，所选点是射线与表面相交的位置。如上所述，它将具有偏差，但是使用小体素，偏差可能足够小。
有一些库可以非常高效地进行凸形状交集计算，其中立方体/椭球体将是其中之一。它们将被高度优化，但我认为无论如何手动执行距离剔除可能都是值得的。您还需要一个库来区分表面交集和一个对象完全在另一个对象内部。
如果您知道您的椭球体与轴对齐，则可以将体素/边缘交集视为一堆2D正方形交叉椭圆问题，要测试的正方形集定义为那些与上层相邻的正方形。这可能会更快。
使这种方法更易于管理的一件事是，您不需要编写所有边缘情况的代码（解决浮点精度问题可能导致交点丢失或重复的工作量很大）。因为这些情况非常罕见，所以它们不会影响您的采样。
甚至可能更快的方法是找到椭圆内的所有体素，然后丢弃所有具有6个邻居的体素...有很多选择。这完全取决于性能的重要性。这比其他建议要慢得多，但如果您想要~1000个点，则表面上大约需要~100,000个体素，因此您可能需要在边界框中使用~1,000,000个体素。但是，即使测试1,000,000个交点在现代计算机上也非常快。
- Bill Sellers
网页内容由stack overflow 提供, 点击上面的可以查看英文原文，
原文链接
原文链接
Ax²+Bxy+Cy²+Dx+Ey+F=0吗？ - Stéphane Laurent