SimpleITK下非单位方向矩阵3D容积CT按欧拉角绕三轴旋转的实现问题

核心原理 #

要实现非单位方向矩阵下的图像自身坐标轴欧拉旋转，核心是区分图像局部坐标系和世界物理坐标系的变换映射：

1. 记图像原方向矩阵为D（3x3，从image.GetDirection()获取后reshape为3x3格式），D_inv为D的逆矩阵，D的作用是将图像局部坐标映射到世界物理坐标
2. 你所需要的XYZ欧拉旋转是图像局部坐标系下的旋转变换，对应旋转矩阵记为R_euler
3. 最终应用到SimpleITK变换的旋转矩阵为R_final = D @ R_euler @ D_inv，通过这个映射可以保证旋转完全沿图像自身的x/y/z轴执行，不受原方向矩阵影响

完整代码实现 #

首先补全缺失的matrix_from_angle辅助函数，再修改原rotation3d函数支持三轴欧拉角输入：

import numpy as np
import SimpleITK as sitk
import matplotlib.pyplot as plt

def matrix_from_angle(axis, theta):
    """
    生成绕指定单位轴旋转指定弧度的旋转矩阵
    :param axis: 旋转轴，0=x轴，1=y轴，2=z轴
    :param theta: 旋转弧度
    :return: 3x3旋转矩阵
    """
    c = np.cos(theta)
    s = np.sin(theta)
    if axis == 0:
        return np.array([[1, 0, 0],
                         [0, c, -s],
                         [0, s, c]])
    elif axis == 1:
        return np.array([[c, 0, s],
                         [0, 1, 0],
                         [-s, 0, c]])
    elif axis == 2:
        return np.array([[c, -s, 0],
                         [s, c, 0],
                         [0, 0, 1]])
    else:
        raise ValueError("axis must be 0, 1, or 2")

# 保留原有已验证的辅助函数
def matrix_from_axis_angle(a):
    """ Compute rotation matrix from axis-angle.
    This is called exponential map or Rodrigues' formula.
    Parameters
    ----------
    a : array-like, shape (4,)
        Axis of rotation and rotation angle: (x, y, z, angle)
    Returns
    -------
    R : array-like, shape (3, 3)
        Rotation matrix
    """
    ux, uy, uz, theta = a
    c = np.cos(theta)
    s = np.sin(theta)
    ci = 1.0 - c
    R = np.array([[ci * ux * ux + c,
                   ci * ux * uy - uz * s,
                   ci * ux * uz + uy * s],
                  [ci * uy * ux + uz * s,
                   ci * uy * uy + c,
                   ci * uy * uz - ux * s],
                  [ci * uz * ux - uy * s,
                   ci * uz * uy + ux * s,
                   ci * uz * uz + c],
                  ])
    return R

def matrix_from_euler_xyz(e):
    """Compute rotation matrix from xyz Euler angles.
    Intrinsic rotations are used to create the transformation matrix
    from three concatenated rotations.
    The xyz convention is usually used in physics and chemistry.
    Parameters
    ----------
    e : array-like, shape (3,)
        Angles for rotation around x-, y'-, and z''-axes (intrinsic rotations,单位：弧度)
    Returns
    -------
    R : array-like, shape (3, 3)
         Rotation matrix
    """
    alpha, beta, gamma = e
    Qx = matrix_from_angle(0, alpha)
    Qy = matrix_from_angle(1, beta)
    Qz = matrix_from_angle(2, gamma)
    R = Qx.dot(Qy).dot(Qz)
    return R

def resample(image, transform):
   """
   基于指定变换对图像进行重采样
   :param image: 输入sitk图像
   :param transform: sitk变换对象
   :return: 变换后的sitk图像
   """
   reference_image = image
   interpolator = sitk.sitkLinear
   default_value = 0
   return sitk.Resample(image, reference_image, transform,
                        interpolator, default_value)

def get_center(img):
   """
   获取3D sitk图像的物理中心坐标
   :param img: 输入sitk图像
   :return: 物理中心坐标元组
   """
   width, height, depth = img.GetSize()
   return img.TransformIndexToPhysicalPoint((int(np.ceil(width/2)),
                                             int(np.ceil(height/2)),
                                             int(np.ceil(depth/2))))

# 修改后的三轴欧拉角旋转函数
def rotation3d(image, theta_x, theta_y, theta_z, show=False):
    """
    沿图像自身x/y/z轴做欧拉角旋转（内旋，顺序x→y'→z''）
    :param image: 输入sitk 3D图像
    :param theta_x: 绕x轴旋转角度（单位：度）
    :param theta_y: 绕y轴旋转角度（单位：度）
    :param theta_z: 绕z轴旋转角度（单位：度）
    :param show: 是否可视化结果切片
    :return: 旋转后的sitk图像
    """
    # 角度转弧度
    alpha = np.deg2rad(theta_x)
    beta = np.deg2rad(theta_y)
    gamma = np.deg2rad(theta_z)
    
    euler_transform = sitk.Euler3DTransform()
    # 设置旋转中心为图像物理中心，避免旋转后偏移
    image_center = get_center(image)
    euler_transform.SetCenter(image_center)

    # 获取原图像方向矩阵，reshape为3x3
    direction = np.array(image.GetDirection()).reshape(3,3)
    # 生成局部坐标系下的欧拉旋转矩阵
    R_euler = matrix_from_euler_xyz([alpha, beta, gamma])
    # 映射到世界坐标系下的变换矩阵
    R_final = direction @ R_euler @ np.linalg.inv(direction)
    # 设置变换矩阵
    euler_transform.SetMatrix(R_final.flatten().tolist())
    # 重采样得到结果
    resampled_image = resample(image, euler_transform)
    
    if show:
        slice_num = int(input("输入要查看的切片索引: "))
        plt.imshow(sitk.GetArrayFromImage(resampled_image)[slice_num], cmap='gray')
        plt.axis('off')
        plt.show()
    return resampled_image

使用示例 #

# 读取你的CT图像
ct_img = sitk.ReadImage("your_ct_path.nii.gz")
# 沿图像自身x轴转5度，y轴转3度，z轴转-2度
rotated_img = rotation3d(ct_img, theta_x=5, theta_y=3, theta_z=-2, show=True)
# 保存结果
sitk.WriteImage(rotated_img, "rotated_ct.nii.gz")

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