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

遗憾的是，这并没有达到预期目的。这就是原因:

定义:

In the z-x-z convention, the x-y-z frame is rotated three times: first about the z-axis by an angle phi; then about the new x-axis by an angle psi; then about the newest z-axis by an angle 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”的第一个参数调整为您的欧拉约定即可。在这种情况下，您获得的结果与我的第一篇文章(使用 scipy.ndimage.rotate)相同。然而我注意到，当前的代码比使用 scipy.ndimage.rotate(40^3 体积为 0.03 秒)快 3 倍(40^3 体积为 0.01 秒)。

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