将欧拉角定义的旋转应用到 3D 图像，python

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

我在这个link之后找到了一个令人满意的解决方案：https://nbviewer.jupyter.org/gist/lhk/f05ee20b5a826e4c8b9bb3e528348688

此方法使用np.meshgrid、scipy.ndimage.map_coordinates。上面的 link 使用了一些第三方库来生成旋转矩阵，但是我使用 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' 的第一个参数适应您的欧拉约定即可。在这种情况下，您获得的结果与我的第一个 post（使用 scipy.ndimage.rotate）相同。但是我注意到当前代码比使用 scipy.ndimage.rotate（40^3 体积为 0.03s）时快 3 倍（40^3 体积为 0.01s）。

希望这对某人有所帮助！

Answer 2

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

Answer 3

我检查代码https://nbviewer.jupyter.org/gist/lhk/f05ee20b5a826e4c8b9bb3e528348688 有一个小错误。测试的图片是正方形的，如果做矩形图片，会遇到一些问题。下面是正确的 2D 和 3D 旋转（请注意，我的示例中使用的欧拉角序列是 'ZYZ'，您应该在使用它之前定义它）：

def rotate_array_2D（数组，方向）：

# create a transformation matrix
angle=orient/180.*np.pi
c=np.cos(angle)
s=np.sin(angle)
mat=np.array([[c,s],[-s,c]])


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

coords = np.meshgrid(ax, ay)

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



# apply transformation
transformed_xy = np.dot(mat, xy)

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


x = x.reshape((dim[1],dim[0]))
y = y.reshape((dim[1],dim[0]))
new_xy = [x,y]

# sample
arrayR = map_coordinates(array, new_xy, order=1).T
return arrayR

def rotate_array_3D（数组，方向）：

rot = orient[0]
tilt = orient[1]
phi = 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('ZYZ', [rot, tilt, phi], 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])) # I test the rotation in 2D and this strange thing can be explained
new_xyz = [x,y,z]

arrayR = map_coordinates(array, new_xyz, order=1).T
return arrayR

将欧拉角定义的旋转应用到 3D 图像，python

Apply rotation defined by Euler angles to 3D image, in python

python

3d

image-processing

euler-angles

image-rotation