如何基于参考平面在指定体素位置切割NIFTI图像？

核心思路 #

你需要先从旋转后的平面点推导出平面方程，再在该平面上生成采样网格，将网格点从世界坐标转回体素坐标，最后通过插值获取体素值，生成切片。同时要修正原代码中NIFTI坐标转换的错误（原代码忽略了affine的平移分量）。

完整实现步骤与代码 #

1. 修正坐标转换逻辑 #

NIFTI的affine是4x4矩阵，完整的体素→世界坐标转换应该用齐次坐标：world_coord = affine @ np.hstack([voxel_coord, 1])，反之世界→体素坐标用affine的逆矩阵：voxel_coord = np.linalg.inv(affine) @ np.hstack([world_coord, 1])。

2. 推导平面方程 #

已知平面上3个点，可计算平面的法向量，进而得到平面方程a(x-x0)+b(y-y0)+c(z-z0)=0，整理为ax + by + cz = d。

3. 生成平面采样网格 #

在平面上生成均匀分布的网格点，确保覆盖足够的图像范围。

4. 插值获取切片数据 #

用scipy.ndimage.map_coordinates实现三线性插值，获取采样点的体素值。

完整代码：

import nibabel as nib
import numpy as np
from scipy.spatial.transform import Rotation
from scipy.ndimage import map_coordinates
import matplotlib.pyplot as plt

def get_custom_slice(image_path, voxel_center, ref_plane_points, rotation_vec):
    # 加载NIFTI数据
    img = nib.load(image_path)
    data = img.get_fdata()
    affine = img.affine
    inv_affine = np.linalg.inv(affine)
    
    # 1. 处理参考平面：旋转后平移到目标体素中心
    # 参考平面点（体素空间）
    ref_plane_voxel = np.array(ref_plane_points)
    # 生成旋转矩阵
    rot_mat = Rotation.from_rotvec(rotation_vec).as_matrix()
    # 旋转参考平面点（体素空间）
    rotated_plane_voxel = (rot_mat @ ref_plane_voxel.T).T
    # 将目标体素中心转成世界坐标
    center_world = affine @ np.hstack([voxel_center, 1])[:3]
    # 旋转后的平面原点（体素空间）转世界坐标，计算平移量
    plane_origin_world = affine @ np.hstack([rotated_plane_voxel[0], 1])[:3]
    translate_vec = center_world - plane_origin_world
    # 平移旋转后的平面点到目标中心（世界空间）
    rotated_translated_plane_world = rotated_plane_voxel @ affine[:3,:3].T + translate_vec
    
    # 2. 推导平面方程：ax + by + cz = d
    p0, p1, p2 = rotated_translated_plane_world[:3]
    # 计算平面法向量
    normal = np.cross(p1 - p0, p2 - p0)
    normal = normal / np.linalg.norm(normal)  # 归一化
    a, b, c = normal
    d = a*p0[0] + b*p0[1] + c*p0[2]
    
    # 3. 生成平面采样网格（世界空间）
    # 定义平面内的两个正交方向向量
    vec1 = p1 - p0
    vec1 = vec1 / np.linalg.norm(vec1)
    vec2 = np.cross(normal, vec1)
    vec2 = vec2 / np.linalg.norm(vec2)
    
    # 设置切片的尺寸（根据图像范围调整）
    slice_size = 64  # 64x64的切片
    step = 2.0  # 采样步长（世界空间单位，比如mm）
    x = np.linspace(-slice_size/2 * step, slice_size/2 * step, slice_size)
    y = np.linspace(-slice_size/2 * step, slice_size/2 * step, slice_size)
    xx, yy = np.meshgrid(x, y)
    
    # 生成网格点的世界坐标
    grid_world = center_world[None, None, :] + xx[..., None]*vec1 + yy[..., None]*vec2
    
    # 4. 世界坐标转回体素坐标
    grid_homogeneous = np.concatenate([grid_world, np.ones(grid_world.shape[:2] + (1,))], axis=-1)
    grid_voxel = (inv_affine @ grid_homogeneous.reshape(-1,4).T).T[:, :3]
    
    # 5. 三线性插值获取切片数据
    slice_data = map_coordinates(data, grid_voxel.T, order=1, mode='constant')
    slice_data = slice_data.reshape(slice_size, slice_size)
    
    return slice_data, rotated_translated_plane_world

# 参数设置
image_path = 'brain_ct.nii'
voxel_center = np.array([-12, -8, 12])  # 目标体素位置
ref_plane_points = np.array([[0,0,0], [1,0,0], [0,1,0], [1,1,0]])  # 原始XY平面参考点
rotation_vec = np.radians([2.4, 0, 0])  # X轴旋转2.4度

# 获取切片
slice_data, plane_points = get_custom_slice(image_path, voxel_center, ref_plane_points, rotation_vec)

# 可视化切片
plt.figure(figsize=(8,8))
plt.imshow(slice_data, cmap='gray')
plt.title('Custom Aligned Slice')
plt.axis('off')
plt.show()

关键细节解释 #

• 坐标转换：必须使用完整的4x4 affine矩阵，包含平移分量，否则会出现位置偏移。
• 平面方程：通过三个点的叉乘计算法向量，确保平面方向正确。
• 采样网格：基于平面内的正交向量生成网格，保证切片的均匀性。
• 插值方式：order=1是三线性插值，平衡精度和速度；mode='constant'处理超出图像范围的点，填充0（可根据需求调整）。

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