如何对Python中坐标矩阵形式的椭球执行欧拉旋转并绘图？

你完全不需要把每个坐标点单独转成向量再处理——我们可以利用NumPy的向量化特性，直接对整个坐标矩阵执行旋转操作，既高效又简洁。下面是具体的实现步骤：

1. 适配原旋转函数到多维坐标 #

你的rot_euler函数目前是针对二维向量数组（形状为(N,3)）设计的，但我们的椭球坐标是三维矩阵（每个坐标矩阵形状为(100,100)）。只需要稍微修改函数，让它支持更高维度的输入，同时保留欧拉旋转的核心逻辑。

2. 修改后的旋转函数（支持多维数组） #

核心思路是用**矩阵乘法运算符@**代替原函数的np.dot，它能自动处理多维数组的广播，批量完成所有点的旋转，无需手动循环每个坐标点：

import numpy as np
from scipy.linalg import expm

def rot_euler_multidim(points, xyz):
    '''
    对多维坐标数组执行欧拉旋转
    参数:
        points: 形状为(3, H, W)的数组，对应x/y/z三个坐标矩阵
        xyz: 欧拉角数组，格式为(roll, pitch, yaw)
    返回:
        旋转后的坐标数组，形状与输入一致
    '''
    for theta, axis in zip(xyz, np.eye(3)):
        # 生成当前轴的旋转矩阵
        R = expm(np.cross(np.eye(3), axis * -theta))
        # 批量旋转所有点：自动广播到所有(H,W)维度的坐标
        points = R @ points
    return points

3. 整合到椭球绘图代码中 #

把旋转逻辑嵌入你的现有代码，只需要几步：

• 将ex、ey、ez堆叠成一个(3, 100, 100)的统一数组
• 调用修改后的旋转函数
• 拆分回旋转后的x/y/z坐标矩阵
• 对比绘制原始椭球和旋转后的椭球

完整代码示例：

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

def rot_euler_multidim(points, xyz):
    for theta, axis in zip(xyz, np.eye(3)):
        R = expm(np.cross(np.eye(3), axis * -theta))
        points = R @ points
    return points

# 初始化3D绘图
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(30, 45)  # 调整视角方便观察旋转效果

# 生成原始椭球坐标
ellipseSteps= 100
coefs = (25, 500, 500)
rx, ry, rz = 1/np.sqrt(coefs)

u = np.linspace(0, 2 * pi, ellipseSteps)
v = np.linspace(0, pi, ellipseSteps)

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))

# 执行欧拉旋转：示例为绕x轴转45°、y轴转30°、z轴转60°
rotation_angles = (pi/4, pi/6, pi/3)
# 堆叠坐标为(3, H, W)的数组
original_points = np.stack([ex, ey, ez], axis=0)
rotated_points = rot_euler_multidim(original_points, rotation_angles)
# 拆分回x/y/z矩阵
ex_rot, ey_rot, ez_rot = rotated_points

# 绘制原始椭球（半透明灰色）和旋转后的椭球（蓝色）
ax.plot_surface(ex, ey, ez, rstride=4, cstride=4, color='gray', alpha=0.3)
ax.plot_surface(ex_rot, ey_rot, ez_rot, rstride=4, cstride=4, color='blue')

plt.show()

为什么这种方式高效？ #

• 完全避免了Python层面的循环遍历每个点，所有计算都在NumPy的C底层执行，速度比逐个处理向量快几个数量级
• 旋转矩阵只需要生成3次（对应三个欧拉轴），而非每个点都生成一次，进一步节省了计算资源

内容的提问来源于stack exchange，提问作者PTTHomps