Эйлерово вращение эллипсоида, выраженное координатными матрицами в питоне

Question

Цель: применить вращение Эйлера к эллипсоиду, затем построить его, используя matplotlib и mplot3d.

Я нашел функцию, которая применяет вращение Эйлера к вектору или массиву векторов:

import numpy as np
from scipy.linalg import expm

def rot_euler(v, xyz):
''' Rotate vector v (or array of vectors) by the euler angles xyz '''
# /4201267/vraschenie-3d-vektora
for theta, axis in zip(xyz, np.eye(3)):
    v = np.dot(np.array(v), expm(np.cross(np.eye(3), axis*-theta)))
return v 

Тем не менее, эллипсоид, который мне нужно повернуть, представлен в виде набора из трех координатных матриц (для построения с использованием ax.plot_surface ()):

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from numpy import pi,sin,cos

fig, ax = plt.subplots(subplot_kw=dict(projection='3d'))
ax.set_aspect('equal','box')
ax.set_xlim3d(-1,1)
ax.set_ylim3d(-1,1)
ax.set_zlim3d(-1,1)
ax.view_init(90,90)

ellipseSteps= 100
diamCoef = 500
widthCoef = 25
coefs = (widthCoef, diamCoef, diamCoef)  # 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 * pi, ellipseSteps)
v = np.linspace(0, pi, ellipseSteps)

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

# Plot:
ax.plot_surface(ex, ey, ez,  rstride=4, cstride=4, color='blue')

plt.show()

Как применить вращение Эйлерак этому объекту?Я думал о том, чтобы просто преобразовать объект из трех координатных матриц в векторы и затем передать его в существующую формулу, но мне приходит в голову, что это может быть вычислительно неэффективно ... что заставляет меня задуматься, можно ли изменить функцию вращения для работыпо координатным матрицам?

Я понимаю, что это, вероятно, довольно тривиальный вопрос, но прошло много лет с тех пор, как я последний раз занимался линейной алгеброй, и буду очень признателен совету эксперта.

Спасибо!

Answer 1

Хотя я не думаю, что вы можете добиться гораздо большего успеха, чем применение поворота по точкам, все же есть место для значительной экономии.

(1) Использование экспоненциальной матрицы для вычисления простых матриц вращения смехотворно расточительно,Намного лучше использовать скалярную экспоненту или синус и косинус

(2) В меньшей степени то же самое относится к использованию перекрестного произведения для тасования.Индексирование здесь предпочтительнее.

(3) Порядок умножения матриц имеет значение.При массовом вращении более трех векторов самое левое умножение должно выполняться последним.

Вместе эти меры ускоряют вычисление в шесть раз:

(вставьте перед последними двумя строками оригиналаскрипт)

from scipy.linalg import block_diag
from timeit import timeit

def rot_euler_better(v, xyz):
    TD = np.multiply.outer(np.exp(1j * np.asanyarray(xyz)), [[1], [1j]]).view(float)
    x, y, z = (block_diag(1, TD[i])[np.ix_(*2*(np.arange(-i, 3-i),))] for i in range(3))
    return v @ (x @ y @ z)

# example
xyz = np.pi * np.array((1/6, -2/3, 1/4))

print("Same result:",
      np.allclose(rot_euler(np.array((*map(np.ravel, (ex, ey, ez)),)).T, xyz),
                  rot_euler_better(np.array((*map(np.ravel, (ex, ey, ez)),)).T, xyz))

print("OP:       ", timeit(lambda: rot_euler(np.array((*map(np.ravel, (ex, ey, ez)),)).T, xyz), number=1000), "ms")
print("optimized:", timeit(lambda: rot_euler_better(np.array((*map(np.ravel, (ex, ey, ez)),)).T, xyz), number=1000), "ms")

ex, ey, ez = map(np.reshape, rot_euler_better(np.array((*map(np.ravel, (ex, ey, ez)),)).T, xyz).T, map(np.shape, (ex, ey, ez)))

Вывод:

---

Same result: True
OP:        2.1019406360574067 ms
optimized: 0.3485010238364339 ms