Примечание, но то, что у вас есть, не является самым общим уравнением для 3-го эллипсоида. Ваше уравнение может быть переписано как
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
, где @ обозначает произведение матрицы на вектор или вектор * .
Вы можете взять свою функцию и построить ееisosurface , но это было бы неоптимальным: для вашей функции вам понадобится сеточное приближение, которое очень дорого сделать с достаточным разрешением, и вам придется выбирать домен для этой выборки с умом.
Вместо этого вы можете выполнить преобразование главной оси для ваших данных, чтобы обобщить параметрический график канонического эллипсоида , который вы сами связали.
Первый шаг - это диагонализация M как M = V @ D @ V.T, где D равно диагональ . Поскольку это реальная симметричная матрица, это всегда возможно, и V является ортогональным . Тогда у нас есть
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).
Итак, вот что мы делаем:
- гарантируют, чтопараметры соответствуют эллипсоиду
- генерация тэты и фи сетки для полярных и азимутальных углов
- вычисление преобразованных координат
[x2, y2, z2] - смещение их назад (на
r2 - r1), чтобы получить [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 масштабирования из этого ответа для сохранения равных осей:
# 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 на функциюсоздать сетку с другой сеткой. Самое интересное, что вы можете взять любые рассеянные точки, лежащие на сфере устройства, и подключить их к определению r2 в функции get_ellipsoid_coordinates (при условии, что этот массив имеет первое измерение размера3), и выходные координаты будут иметь ту же форму, но они будут преобразованы в фактический эллипсоид.
Вы также можете использовать другие библиотеки для визуализации поверхности, например, Mayavi, где вы можете либо нарисовать поверхностьмы только что вычислили или сравнили его с изоповерхностью , которая там встроена.