使用Matplotlib绘制椭圆体

Question

有没有人有绘制椭圆体的示例代码？matplotlib现场有一个球体,但椭圆体没有.我想策划

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

其中c是一个定义椭球的常数(如10).我尝试了这meshgrid(x,y)条路线,重新设计了等式,所以z就在一边,但这sqrt是一个问题.该matplotlib球体例子使用角度u,v,但我不知道如何工作,对于椭球.

Answer 1

以下是通过球坐标来完成的方法:

from mpl_toolkits.mplot3d import Axes3D
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的画廊.

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为中心并且想要绘制它.

Answer 3

如果您有一个由任意协方差矩阵cov和 offset指定的椭球体，则可以通过向量化操作来执行@minillinim 答案bias的简单版本。

从...开始

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)，对于广播运营商而言，这使其成为 3x1 向量的 100x100 数组@。s @ sphere大小相同，所以挤出最后一个单位轴，适合用 进行加法广播bias。最后，ellipsoid.transpose(2, 0, 1)有 shape (3, 100, 100)，它可以星形展开为三个独立的 x、y 和 z 值数组，并调用plot_surface。