从MRI DICOM生成3D仿射变换矩阵：RTDose重采样失败求助

AI 小北

2026-6-22

问题描述

我有DICOM格式的RTDose与MRI数据，二者维度不同、实际间距也不兼容。参考DICOM仿射定义文档尝试构建3D仿射变换矩阵，将RTDose重采样至MRI的间距，但可视化切片后发现变换未生效。实现代码如下：

#read mri
mr_file_paths = [os.path.join(path_MR, f) for f in os.listdir(path_MR) if f.endswith('.dcm')]        
mr_data = [pydicom.dcmread(p) for p in mr_file_paths] 
mr_dtype = mr_data[0].pixel_array.dtype
mr = np.array([p.pixel_array for p in mr_data], dtype=mr_dtype)

#read rtdose
rtd_file_path  = [os.path.join(path_RTD, f) for f in os.listdir(path_RTD) if f.endswith('.dcm')][0]
rtd_data = pydicom.dcmread(rtd_file_path)
rtd = np.float64(rtd_data.pixel_array)

#get patient's position, orientation and mri pixel spacing
mr_ippf, mr_ippl = np.array(mr_data[0].ImagePositionPatient), np.array(mr_data[-1].ImagePositionPatient)
mr_iop = np.array(mr_data[0].ImageOrientationPatient)
mr_spacing = np.array(mr_data[0].PixelSpacing)

Δr, Δc = mr_spacing[0], mr_spacing[1]
F11, F21, F31 = mr_iop[3:6]
F12, F22, F32 = mr_iop[0:3]

k1, k2, k3 = (mr_ippl - mr_ippf) / (len(mr)-1)
Dx, Dy, Dz = mr_ippf
#build affine transformation matrix
affine = np.array([
    [F11*Δr,    F12*Δc, k1, Dx  ],
    [F21*Δr,    F22*Δc, k2, Dy  ],
    [F31*Δr,    F32*Δc, k3, Dz  ],
    [0,         0,      0,  1   ],
], dtype=np.float64)
#apply affine transformation
rtd = ndimage.affine_transform(rtd, affine, output_shape=mr.shape)

解决建议

纠正变换矩阵的方向：scipy.ndimage.affine_transform需要的是从输出空间到输入空间的逆变换矩阵，你现在传入的是MRI体素到世界坐标的正向矩阵，完全逻辑颠倒。正确做法是：
1. 分别计算RTDose和MRI的世界坐标仿射矩阵（体素坐标→世界坐标）
2. 计算变换矩阵：transform = np.linalg.inv(affine_mr) @ affine_rtd，该矩阵实现RTDose体素坐标→世界坐标→MRI体素坐标的映射，这才是affine_transform所需的参数。
核对维度顺序一致性：RTDose的pixel_array维度通常为(Z, Y, X)，要确认和你堆叠得到的MRI数组mr.shape（(Z, 高度, 宽度)）的维度顺序完全匹配，维度错位会直接导致变换失效。
避免手动构建仿射矩阵：手动计算容易出现IOP参数顺序、层间距计算等错误，建议用成熟工具生成标准仿射：
- MRI序列：用nibabel.dicom.read_dicom_dir(path_MR)读取后直接获取.affine属性
- RTDose：用nibabel.dicom.read_dicom_file(rtd_file_path)读取后获取.affine属性
验证世界坐标范围：变换前打印两者的世界坐标边界（通过仿射矩阵计算体素原点和对角点的世界坐标），确认数据处于同一患者坐标系内，排除DICOM文件定位参数错误的可能。
优化插值与边缘处理：剂量数据重采样推荐使用三次样条插值（order=3），同时设置边缘填充模式避免异常值，示例代码：
```
rtd_resampled = ndimage.affine_transform(rtd, transform, output_shape=mr.shape, order=3, mode='constant')
```

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