python中坐标矩阵表示的椭球欧拉旋转

Question

Objective: 对椭圆体应用欧拉旋转，然后使用 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 '''
# 
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

python中坐标矩阵表示的椭球欧拉旋转

Euler rotation of ellipsoid expressed by coordinate matrices in python

python

numpy

matplotlib

rotation