使用Matplotlib绘制椭球形图形

Question

有没有人有绘制椭球体的示例代码？在matplotlib网站上有一个绘制球体的示例代码，但没有绘制椭球体的。我正在尝试绘制。

x**2 + 2*y**2 + 2*z**2 = c

其中c是一个常数（比如10），用于定义椭球体形状。我尝试了使用meshgrid(x,y)，重构了方程使得z在一侧，但是sqrt成为一个问题。 matplotlib中的球体示例使用角度u,v，但我不确定如何将其应用于椭球体。

- BBSysDyn

---

@tillsten：您提供的参考文献描述了如何绘制一个椭圆形，而问题是关于椭球体的，它是椭圆的三维等价物。 - Eric O. Lebigot

EOL: 你说得对，我会删除我的评论。 - tillsten

Answer 1

以下是如何通过球坐标系完成此操作：

# from mpl_toolkits.mplot3d import Axes3D  # Not needed with Matplotlib 3.6.3
import matplotlib.pyplot as plt
import numpy as np

fig = plt.figure(figsize=plt.figaspect(1))  # Square figure
ax = fig.add_subplot(111, projection='3d')

coefs = (1, 2, 2)  # Coefficients in a0/c x**2 + a1/c y**2 + a2/c z**2 = 1 
# Radii corresponding to the coefficients:
rx, ry, rz = 1/np.sqrt(coefs)

# Set of all spherical angles:
u = np.linspace(0, 2 * np.pi, 100)
v = np.linspace(0, np.pi, 100)

# Cartesian coordinates that correspond to the spherical angles:
# (this is the equation of an ellipsoid):
x = rx * np.outer(np.cos(u), np.sin(v))
y = ry * np.outer(np.sin(u), np.sin(v))
z = rz * np.outer(np.ones_like(u), np.cos(v))

# Plot:
ax.plot_surface(x, y, z,  rstride=4, cstride=4, color='b')

# Adjustment of the axes, so that they all have the same span:
max_radius = max(rx, ry, rz)
for axis in 'xyz':
    getattr(ax, 'set_{}lim'.format(axis))((-max_radius, max_radius))

plt.show()

生成的图形与...类似。

上述程序实际上生成了一个外观更好的“正方形”图形。

这个解决方案受到Matplotlib画廊中示例的强烈启发。

- Eric O Lebigot

---

没错，那就是问题所在。我升级到了1.1.0版本，现在一切都正常了。 - BBSysDyn

rx, ry, rz = 1/np.sqrt(coefs) - derchambers

如果我的半长轴为 a 和 c，半短轴为 b，那么系数是不是应该是 (1/a2, 1/b2, 1/c**2)？ - Alex Punnen

当然可以直接使用 rx, ry, rz = a, b, c 进行赋值（不需要系数 coefs）。 - Eric O. Lebigot

对于像我一样最初不理解np.outer如何工作的人，可以使用以下代码替代：u = np.linspace(0, 2 * np.pi, 100); v = np.linspace(0, np.pi, 100); u,v = np.meshgrid(u,v,indexing='ij'); x = rx * np.cos(u) * np.sin(v); y = ry * np.sin(u) * np.sin(v); z = rz * np.cos(v)。 - mhdadk

Answer 2

在 EOL 的回答基础上进行拓展。有时你会遇到以矩阵格式表示的椭球体：

A 和 c，其中 A 是椭球体的矩阵，而 c 则是代表椭球体中心的向量。

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

# your ellispsoid and center in matrix form
A = np.array([[1,0,0],[0,2,0],[0,0,2]])
center = [0,0,0]

# find the rotation matrix and radii of the axes
U, s, rotation = linalg.svd(A)
radii = 1.0/np.sqrt(s)

# now carry on with EOL's answer
u = np.linspace(0.0, 2.0 * np.pi, 100)
v = np.linspace(0.0, np.pi, 100)
x = radii[0] * np.outer(np.cos(u), np.sin(v))
y = radii[1] * np.outer(np.sin(u), np.sin(v))
z = radii[2] * np.outer(np.ones_like(u), np.cos(v))
for i in range(len(x)):
    for j in range(len(x)):
        [x[i,j],y[i,j],z[i,j]] = np.dot([x[i,j],y[i,j],z[i,j]], rotation) + center

# plot
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_wireframe(x, y, z,  rstride=4, cstride=4, color='b', alpha=0.2)
plt.show()
plt.close(fig)
del fig

所以，这里没有太多新东西，但如果您有一个矩阵形式的椭球体，它被旋转并且可能不位于 0,0,0 的中心位置，并且想要绘制它，那么这个信息是有帮助的。

- minillinim

---

你好，为什么需要计算奇异值的倒数平方根来获取半径，而不是直接使用平方根呢？我的意思是，为什么半径等于1.0/np.sqrt(s)，而不是np.sqrt(s)。我认为后者可以得到正确的半径。 - isti_spl

你根本不需要这些复杂的数学，尤其是奇异值分解。协方差矩阵的整个意义在于，如果你用它乘以一个球体，你会得到一个椭圆。 - Mad Physicist

我已经发布了一个回答，演示了我的意思。 - Mad Physicist

Answer 3

如果你有一个由任意协方差矩阵cov和偏移bias指定的椭球体，你可以通过向量化操作来执行@minillinim's answer的简化版本。

从以下开始：

u = np.linspace(0, 2 * np.pi, 100)
v = np.linspace(0, np.pi, 100)

制作一个单位球

x = np.outer(np.cos(u), np.sin(v))
y = np.outer(np.sin(u), np.sin(v))
z = np.outer(np.ones_like(u), np.cos(v))
sphere = np.stack((x, y, z), axis=-1)[..., None]

计算标准差矩阵，该矩阵具有与协方差相同的旋转，但通过特征值的平方根进行缩放。

e, v = np.linalg.eig(cov)
s = v @ np.diag(np.sqrt(e)) @ v.T

现在将球体进行变换：

ellipsoid = (s @ sphere).squeeze(-1) + bias

你可以像之前一样漂亮地绘制结果。

ax.plot_surface(*ellipsoid.transpose(2, 0, 1), rstride=4, cstride=4, color='b', alpha=0.75)

供参考，u，v 的形状为(100,)，这使得x，y，z 成为(100, 100)的数组。sphere 是(100, 100, 3, 1)，就广播运算符@而言，它是一个100x100的3x1向量数组。s @ sphere具有相同的大小，因此挤压掉最后一个单位轴使其适合与bias进行广播加法。最后，ellipsoid.transpose(2, 0, 1)的形状为(3, 100, 100)，可以将其作为x、y和z值的三个单独数组扩展到plot_surface的调用中。

- Mad Physicist

---

@Translunar。矩阵上的任何标量乘法都是有效的。 - Mad Physicist

然而，在这种情况下，你可以得到它，因为我们的“数据集”是球面上的点。当你在上面执行R'时，你会得到相同的点集，所以你的变换C = RSS'R'实际上是变换RSS'。这几乎就像应该是的RS，只是你将数据两次缩放S而不是一次。这就是为什么Translunar抱怨协方差椭球没有适当地按比例缩放到数据的原因。因此，这种方法确实给出了一个几乎具有正确形状的椭球，但正如我所说，这是一个小技巧。 - Jagerber48

@Jagerber48。你是完全正确的。我修正了答案以纠正错误。 - Mad Physicist

@Jagerber48。你说得完全正确。我已经修正了答案，以纠正这个错误。 - undefined

@Translunar。是的。虽然晚了几年，但现在更新后的答案应该是正确的。 - Mad Physicist