如何基于DICOM可形变配准文件对CT扫描数据执行形变处理

核心逻辑说明 #

你已经提取到的dx、dy、dz是形变矢量场，本质是目标空间每个坐标点对应源CT空间的坐标偏移量，形变操作本质就是通过矢量场计算目标空间每个点在源CT上的对应位置，再通过插值得到形变后的CT值。

前置依赖 #

需要额外安装scipy用于3D插值：pip install scipy

完整实现代码 #

import numpy as np
from scipy.ndimage import map_coordinates

# ---------------------- 前置准备 ----------------------
# 1. 读取待形变的CT数据，确保其形状为 (nz, ny, nx) 和形变场维度完全对齐
# ct_arr = 你从DICOM序列读取得到的CT 3D numpy数组
# 2. 偏移单位转换：如果dx/dy/dz是物理空间偏移（单位mm），需要先转成体素偏移
# 可从CT DICOM的PixelSpacing、SliceThickness字段读取体素间距：
# spacing_x, spacing_y = ct_ds.PixelSpacing
# spacing_z = ct_ds.SliceThickness
# dx_vox = dx / spacing_x
# dy_vox = dy / spacing_y
# dz_vox = dz / spacing_z
# 如果dx/dy/dz已经是体素偏移，直接赋值即可
dx_vox = dx
dy_vox = dy
dz_vox = dz

# ---------------------- 生成采样坐标 ----------------------
# 生成目标空间的原生坐标网格，indexing='ij' 适配 (nz, ny, nx) 的维度顺序
z_grid, y_grid, x_grid = np.meshgrid(np.arange(nz), np.arange(ny), np.arange(nx), indexing='ij')
# 叠加偏移量，得到源CT上的采样坐标
sample_z = z_grid + dz_vox
sample_y = y_grid + dy_vox
sample_x = x_grid + dx_vox

# ---------------------- 插值得到形变后CT ----------------------
# 将采样坐标整理成map_coordinates要求的格式：(3, 采样点总数)
coords = np.vstack([sample_z.ravel(), sample_y.ravel(), sample_x.ravel()])
# 线性插值，边界外的值填充为空气CT值-1000
# 如果是形变分割标签，把order改为0使用最近邻插值，避免生成非整数标签值
warped_ct = map_coordinates(ct_arr, coords, order=1, mode='constant', cval=-1000)
# 恢复成和原CT相同的形状
warped_ct = warped_ct.reshape(nz, ny, nx)

注意事项 #

• 坐标顺序校验：如果运行后形变结果完全错位，优先检查网格的索引方式，部分DICOM形变场的维度顺序可能为(nx, ny, nz)，需要对应调整reshape和meshgrid的参数。
• 形变方向校验：如果形变后结果和预期相反，说明你的形变场是固定图像到浮动图像的映射，此时如果要把浮动图像形变到固定图像空间，需要先计算形变场的逆变换，或调整偏移量的正负号。

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