在python中将由欧拉角定义的旋转应用于3D图像答案

【问题标题】：Apply rotation defined by Euler angles to 3D image, in python在python中将由欧拉角定义的旋转应用于3D图像
【发布时间】：2020-05-01 10:36:21
【问题描述】：

我正在处理 3D 图像，并且必须根据“zxz”约定中的欧拉角（phi、psi、theta）旋转它们（这些欧拉角是数据集的一部分，因此我必须使用该约定）。我发现函数 scipy.ndimage.rotate 在这方面似乎很有用。

arrayR = scipy.ndimage.rotate(array , phi, axes=(0,1), reshape=False)
arrayR = scipy.ndimage.rotate(arrayR, psi, axes=(1,2), reshape=False)
arrayR = scipy.ndimage.rotate(arrayR, the, axes=(0,1), reshape=False)

遗憾的是，这并没有达到预期的效果。这就是为什么：

定义：

在 z-x-z 约定中，x-y-z 框架旋转了 3 次：首先绕 z 轴成角度 phi；然后关于 new x 轴角度 psi;然后关于最新 z 轴的角度 theta。

但是，使用上面的代码，旋转总是相对于原始轴。这就是为什么获得的旋转不正确的原因。如定义中所述，有人建议获得正确的旋转吗？

换句话说，在当前的“zxz”约定中，旋转是固有的（围绕旋转坐标系 XYZ 的轴旋转，与移动体一体，在每次基本旋转后都会改变其方向）。如果我使用上面的代码，旋转是外在的（围绕原始坐标系的轴 xyz 旋转，假设它保持不动）。我需要一种在 python 中进行外部旋转的方法。

【问题讨论】：

标签： python image-processing 3d image-rotation euler-angles

【解决方案1】：

我通过这个链接找到了一个令人满意的解决方案：https://nbviewer.jupyter.org/gist/lhk/f05ee20b5a826e4c8b9bb3e528348688

此方法使用 np.meshgrid、scipy.ndimage.map_coordinates。上面的链接使用一些第三方库来生成旋转矩阵，但是我使用 scipy.spatial.transform.Rotation。此函数允许定义内在和外在旋转：参见 scipy.spatial.transform.Rotation.from_euler 的描述。

这是我的功能：

import numpy as np
from scipy.spatial.transform import Rotation as R
from scipy.ndimage import map_coordinates

# Rotates 3D image around image center
# INPUTS
#   array: 3D numpy array
#   orient: list of Euler angles (phi,psi,the)
# OUTPUT
#   arrayR: rotated 3D numpy array
# by E. Moebel, 2020
def rotate_array(array, orient):
    phi = orient[0]
    psi = orient[1]
    the = orient[2]

    # create meshgrid
    dim = array.shape
    ax = np.arange(dim[0])
    ay = np.arange(dim[1])
    az = np.arange(dim[2])
    coords = np.meshgrid(ax, ay, az)

    # stack the meshgrid to position vectors, center them around 0 by substracting dim/2
    xyz = np.vstack([coords[0].reshape(-1) - float(dim[0]) / 2,  # x coordinate, centered
                     coords[1].reshape(-1) - float(dim[1]) / 2,  # y coordinate, centered
                     coords[2].reshape(-1) - float(dim[2]) / 2])  # z coordinate, centered

    # create transformation matrix
    r = R.from_euler('zxz', [phi, psi, the], degrees=True)
    mat = r.as_matrix()

    # apply transformation
    transformed_xyz = np.dot(mat, xyz)

    # extract coordinates
    x = transformed_xyz[0, :] + float(dim[0]) / 2
    y = transformed_xyz[1, :] + float(dim[1]) / 2
    z = transformed_xyz[2, :] + float(dim[2]) / 2

    x = x.reshape((dim[1],dim[0],dim[2]))
    y = y.reshape((dim[1],dim[0],dim[2]))
    z = z.reshape((dim[1],dim[0],dim[2])) # reason for strange ordering: see next line

    # the coordinate system seems to be strange, it has to be ordered like this
    new_xyz = [y, x, z]

    # sample
    arrayR = map_coordinates(array, new_xyz, order=1)

注意：您也可以将此函数用于内部旋转，只需将“from_euler”的第一个参数调整为您的欧拉约定。在这种情况下，您获得的结果与我的第一篇文章（使用 scipy.ndimage.rotate）相同。但是我注意到当前的代码比使用 scipy.ndimage.rotate 时快 3 倍（40^3 体积为 0.01 秒）（40^3 体积为 0.03 秒）。

希望这会对某人有所帮助！

【讨论】：

【解决方案2】：

您的first帖子中的“axes”参数似乎有些混乱。要绕 x 轴旋转，旋转平面将是 yz 平面，这意味着您的“轴”参数应设置为 (1,2)。第一次和第三次旋转大概是围绕 x 和 z 轴进行的。但是，您的两个旋转都在 xy 平面上。这些可能是您答案差异背后的原因吗？我不相信您对新轴和原始轴的解释。对“旋转”函数的独立调用无法访问任何形式或形状的旧数据。它只看到新的轴和旋转的数组。

【讨论】：