将Euler角定义的旋转应用于python中的3D图像

Question

我正在处理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框架旋转三圈：第一绕z轴成角度phi；然后关于new x轴角度psi;然后绕最新 z轴旋转一个角度theta。

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

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

Answer 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”的第一个参数修改为您的Euler约定即可。在这种情况下，您获得的效果与我的第一篇文章相同（使用scipy.ndimage.rotate）。但是我注意到，当前代码比使用scipy.ndimage.rotate（40 ^ 3体积为0.03s）快3倍（对于40 ^ 3体积为0.01s）。

希望这会帮助某人！