我如何在由二次方程定义的椭球的matplotlib中制作3D图?

问题描述 投票:1回答:1

我有一个椭球的通式:

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.
python matplotlib 3d ellipse quadratic-curve
1个回答
0
投票

您拥有的不是3D椭球的最通用方程。

假设我们得到一个总分r0 = [x0, y0, z0],则有

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

where

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.TD。由于它是一个实对称矩阵,因此这始终是可能的,并且diagonalV。然后我们有

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 d1d2d3都严格为正。如果可以保证,则标准形式的半长轴为c2sqrt(c2/d1)sqrt(c2/d2)

所以这就是我们的工作:

  1. 确保参数对应于椭球
  2. 为极角和方位角生成theta和phi网格
  3. 计算变换后的坐标sqrt(c2/d3)
  4. 将它们移回(按[x2, y2, z2])以获得r2 - r1
  5. 将坐标转换为[x1, y1, z1]以获得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()

您可以将函数的figure of an ellipsoid that's flattened and rotated around the z axisn_theta关键字参数更改为生成具有不同网格的网格。有趣的是,您可以获取单位球上的any个分散点,并将其插入到函数n_phi中的r2定义中(只要此数组的初次尺寸为3)。 ,并且输出坐标将具有相同的形状,但是会将它们转换为实际的椭球。

[您还可以使用其他库来可视化表面,例如mayavi,您可以在其中绘制我们刚刚计算出的表面,也可以在其中内置get_ellipsoid_coordinates

© www.soinside.com 2019 - 2024. All rights reserved.