如何在matplotlib中绘制由二次方程定义的椭球体的三维图？

Question

我有一个椭球体的一般公式：

A*x**2 + C*y**2 + D*x + E*y + B*x*y + F + G*z**2 = 0

其中A、B、C、D、E、F、G是常数因子。你知道吗

如何在matplotlib中将此公式绘制为3D绘图？（最好是线框。）

我看到了example，但它是参数形式的，我不知道如何在代码中放置z坐标。有没有一种方法可以保持通用的形式，而不使用参数形式来绘制它？你知道吗

我开始把它放进这样的代码里：

from mpl_toolkits import mplot3d
%matplotlib notebook
import numpy as np
import matplotlib.pyplot as plt

def f(x, y):
    return ((A*x**2 + C*y**2 + D*x + E*y + B*x*y + F))

def f(z):
    return G*z**2

x = np.linspace(-2200, 1850, 30)
y = np.linspace(-100, 60, 30)
z = np.linspace(-100, 60, 30)

X, Y, Z = np.meshgrid(x, y, z)

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.plot_wireframe(X, Y, Z, rstride=10, cstride=10)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z');

我有个错误：

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-1-95b1296ae6a4> in <module>()
     18 fig = plt.figure()
     19 ax = fig.add_subplot(111, projection='3d')
---> 20 ax.plot_wireframe(X, Y, Z, rstride=10, cstride=10)
     21 ax.set_xlabel('x')
     22 ax.set_ylabel('y')

C:\Program Files (x86)\Microsoft Visual Studio\Shared\Anaconda3_64\lib\site-packages\mpl_toolkits\mplot3d\axes3d.py in plot_wireframe(self, X, Y, Z, *args, **kwargs)
   1847         had_data = self.has_data()
   1848         if Z.ndim != 2:
-> 1849             raise ValueError("Argument Z must be 2-dimensional.")
   1850         # FIXME: Support masked arrays
   1851         X, Y, Z = np.broadcast_arrays(X, Y, Z)

ValueError: Argument Z must be 2-dimensional.

Answer 1

旁注，但你所得到的不是一个三维椭球体最一般的方程。你的方程式可以改写为

A*x**2 + C*y**2 + D*x + E*y + B*x*y = - G*z**2 - F,

这意味着对于z的每个值，实际上得到了一个不同级别的2d椭圆，切片相对于z = 0平面是对称的。这表明你的椭球不是一般的，它有助于检查结果，以确保我们得到的是有意义的。你知道吗

假设我们取一个一般点r0 = [x0, y0, z0]，那么

r0 @ M @ r0 + b0 @ r0 + c0 == 0

在哪里

M = [ A    B/2    0
     B/2    C     0
      0     0     G],
b0 = [D, E, 0],
c0 = F

其中^{} stands for matrix-vector or vector-vector product。你知道吗

你可以使用你的函数和plot its isosurface，但那是次优的：你需要一个网格化的函数近似值，这对于足够的分辨率是非常昂贵的，而且你必须明智地选择这个采样的域。你知道吗

相反，您可以对数据执行principal axis transformation，以泛化您自己链接的parametric plot of a canonical ellipsoid。你知道吗

第一步是将M对角化为M = V @ D @ V.T，其中D是diagonal。因为它是实对称矩阵，所以这总是可能的，而且V就是orthogonal。那我们有

r0 @ V @ D @ V.T @ r0 + b0 @ r0 + c0 == 0

我们可以重组为

(V.T @ r0) @ D @ (V.T @ r0) + b0 @ V @ (V.T @ r0) + c0 == 0

它激发了辅助坐标r1 = V.T @ r0和向量b1 = b0 @ V的定义，我们得到

r1 @ D @ r1 + b1 @ r1 + c0 == 0.

因为D是一个对称矩阵，特征值d1, d2, d3在它的对角线上，上面是方程

d1 * x1**2 + d2 * x2**2 + d3 * x3**3 + b11 * x1 + b12 * x2 + b13 * x3 + c0 == 0

其中r1 = [x1, x2, x3]和b1 = [b11, b12, b13]。你知道吗

剩下的就是从r1切换到r2，这样我们就去掉了线性项：

d1 * (x1 + b11/(2*d1))**2 + d2 * (x2 + b12/(2*d2))**2 + d3 * (x3 + b13/(2*d3))**2 - b11**2/(4*d1) - b12**2/(4*d2) - b13**2/(4*d3) + c0 == 0

所以我们定义

r2 = [x2, y2, z2]
x2 = x1 + b11/(2*d1)
y2 = y1 + b12/(2*d2)
z2 = z1 + b13/(2*d3)
c2 = b11**2/(4*d1) b12**2/(4*d2) b13**2/(4*d3) - c0.

为了这些我们终于有了

d1 * x2**2 + d2 * y2**2 + d3 * z2**2 == c2,
d1/c2 * x2**2 + d2/c2 * y2**2 + d3/c2 * z2**2 == 1

这是二阶曲面的标准形式。为了使它有意义地对应于椭球体，我们必须确保d1、d2、d3和c2都是严格正的。如果保证了这一点，则标准形式的半长轴是sqrt(c2/d1)、sqrt(c2/d2)和sqrt(c2/d3)。你知道吗

我们要做的是：

1. 确保参数对应于椭球体
2. 生成极角和方位角的θ和φ网格
3. 计算变换后的坐标[x2, y2, z2]
4. 把它们移回去（按r2 - r1）得到[x1, y1, z1]
5. 将坐标转换回V得到r0，即我们感兴趣的实际[x, y, z]坐标。你知道吗

下面是我如何实现这一点：

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

def get_transforms(A, B, C, D, E, F, G): 
    """ Get transformation matrix and shift for a 3d ellipsoid 

    Assume A*x**2 + C*y**2 + D*x + E*y + B*x*y + F + G*z**2 = 0, 
    use principal axis transformation and verify that the inputs 
    correspond to an ellipsoid. 

    Returns: (d, V, s) tuple of arrays 
        d: shape (3,) of semi-major axes in the canonical form 
           (X/d1)**2 + (Y/d2)**2 + (Z/d3)**2 = 1 
        V: shape (3,3) of the eigensystem 
        s: shape (3,) shift from the linear terms 
    """ 

    # construct original matrix 
    M = np.array([[A, B/2, 0], 
                  [B/2, C, 0], 
                  [0, 0, G]]) 
    # construct original linear coefficient vector 
    b0 = np.array([D, E, 0]) 
    # constant term 
    c0 = F 

    # compute eigensystem 
    D, V = np.linalg.eig(M) 
    if (D <= 0).any(): 
        raise ValueError("Parameter matrix is not positive definite!") 

    # transform the shift 
    b1 = b0 @ V 

    # compute the final shift vector 
    s = b1 / (2 * D) 

    # compute the final constant term, also has to be positive 
    c2 = (b1**2 / (4 * D)).sum() - c0 
    if c2 <= 0: 
        print(b1, D, c0, c2) 
        raise ValueError("Constant in the canonical form is not positive!")

    # compute the semi-major axes 
    d = np.sqrt(c2 / D) 

    return d, V, s 

def get_ellipsoid_coordinates(A, B, C, D, E, F, G, n_theta=20, n_phi=40): 
    """Compute coordinates of an ellipsoid on an ellipsoidal grid 

    Returns: x, y, z arrays of shape (n_theta, n_phi) 
    """ 

    # get canonical grid 
    theta,phi = np.mgrid[0:np.pi:n_theta*1j, 0:2*np.pi:n_phi*1j] 
    r2 = np.array([np.sin(theta) * np.cos(phi), 
                   np.sin(theta) * np.sin(phi), 
                   np.cos(theta)]) # shape (3, n_theta, n_phi) 

    # get transformation data 
    d, V, s = get_transforms(A, B, C, D, E, F, G)  # could be *args I guess 

    # shift and transform back the coordinates 
    r1 = d[:,None,None]*r2 - s[:,None,None]  # broadcast along first of three axes
    r0 = (V @ r1.reshape(3, -1)).reshape(r1.shape)  # shape (3, n_theta, n_phi) 

    return r0  # unpackable to x, y, z of shape (n_theta, n_phi)

下面是一个椭球的例子，证明了它的有效性：

A,B,C,D,E,F,G = args = 2, -1, 2, 3, -4, -3, 4 
x,y,z = get_ellipsoid_coordinates(*args) 
print(np.allclose(A*x**2 + C*y**2 + D*x + E*y + B*x*y + F + G*z**2, 0))  # True

从这里开始的实际策划是微不足道的。使用3d scaling hack from this answer保持轴相等：

# create 3d axes
fig = plt.figure() 
ax = fig.add_subplot(111, projection='3d')

# plot the data
ax.plot_wireframe(x, y, z)
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')

# scaling hack
bbox_min = np.min([x, y, z])
bbox_max = np.max([x, y, z])
ax.auto_scale_xyz([bbox_min, bbox_max], [bbox_min, bbox_max], [bbox_min, bbox_max])

plt.show()

结果如下：

把它旋转一圈就可以很好地看到，表面确实是反射对称的，相对于z = 0平面，这从方程中可以明显看出。你知道吗

您可以更改函数的n_theta和n_phi关键字参数，以生成具有不同网格的网格。有趣的是，你可以把单位球上的任何散乱点插入到函数get_ellipsoid_coordinates中的r2的定义中（只要这个数组的第一个维度是3），输出的坐标将具有相同的形状，但它们将被转换成实际的椭球体。你知道吗

您还可以使用其他库来可视化曲面，例如mayavi，您可以在其中绘制我们刚刚计算的曲面，或者compare it with an isosurface，它内置在那里。你知道吗