4. 图像几何与空间信息
在 SimpleITK 中,一张医学影像不只是像素矩阵,它同时携带一套完整的空间几何描述。像素值回答“这里有多亮”,几何信息回答“这个像素在患者或世界坐标系的哪个位置”。如果忽略几何信息,配准、重采样、ROI 提取、多模态融合都会出现看似正确、实则错位的严重问题。本节围绕四个核心量——Origin、Spacing、Direction,以及索引到物理坐标的转换函数 TransformIndexToPhysicalPoint——说明它们如何共同定义影像的空间位置,并给出坐标系对齐的实用方法。
四个几何量的定义
SimpleITK 使用“索引坐标 → 物理坐标”的仿射映射描述图像几何。设索引为 (i, j, k),则物理点 p 满足:
p = Origin + Direction · (Spacing ⊙ (i, j, k))
其中 ⊙ 表示逐元素相乘,Direction 是 3×3 方向余弦矩阵。
Origin 是索引 (0,0,0) 对应的物理坐标,单位通常为毫米。它并不一定是图像“左上角”在解剖上的起点,而是该图像自身坐标系的原点。Spacing 是三个轴向的体素间距,即相邻索引在物理空间中的距离,例如 CT 常见的 (0.7, 0.7, 1.0) mm。Direction 是 3×3 矩阵,每一列表示一个索引轴在物理坐标系中的方向单位向量。当 Direction 为单位矩阵时,索引轴与物理轴平行;当图像被旋转采集或经过重采样后,Direction 往往不再是单位阵。
TransformIndexToPhysicalPoint 的作用
TransformIndexToPhysicalPoint 把连续索引坐标映射为物理坐标,是几何计算的入口。它接受浮点索引,因此可以表示体素之间的亚像素位置。例如要取第 100 行、第 80 列、第 20 层的体素中心物理坐标:
import SimpleITK as sitk
img = sitk.ReadImage("ct.nii.gz")
idx = (100, 80, 20)
p = img.TransformIndexToPhysicalPoint(idx)
print(p) # 例如 (-123.4, -87.2, 45.6)
# 反向:物理坐标 → 连续索引
i = img.TransformPhysicalPointToContinuousIndex(p)
print(i) # 接近 (100.0, 80.0, 20.0)
这两个函数互为逆映射。实际工程中,标注点、探针位置、配准变换后的目标点通常以物理坐标存储,而访问像素必须回到索引,因此成对使用非常常见。需要注意:TransformIndexToPhysicalPoint 使用的是体素中心约定,索引 (0,0,0) 对应第一个体素的中心,而不是它的角点。
坐标系对齐的常见场景
坐标系对齐的本质,是让不同图像或不同来源的点落在同一物理参考下。典型场景有三类。
- 多模态同机扫描:CT 与 PET 若来自同一设备同一检查,Origin、Spacing、Direction 可能一致,此时索引可直接对应;但仍应显式比较几何量,不能假设。
- 不同分辨率图像:高分辨率 CT 与低分辨率 PET 的 Spacing 不同,同一物理点对应不同索引,必须通过物理坐标转换后再采样。
- 经过配准或重采样:变换会改变 Origin 与 Direction,配准后图像通常以参考图像几何为准,需要重新读取几何量再计算。
下面这张图展示了索引空间与物理空间之间的双向映射关系,以及 Origin、Spacing、Direction 在映射中的位置。图中同时标出一个示例体素从索引到物理坐标的转换路径,便于对照理解。
坐标系对齐的实操方法
当需要把一张图像对齐到另一张参考图像时,最稳妥的做法是显式构造目标几何,而不是手动改 Origin。SimpleITK 提供了 ResampleImageFilter,配合参考图像几何即可完成对齐。下面示例把 moving 图像重采样到 fixed 图像的几何网格上:
fixed = sitk.ReadImage("ct.nii.gz")
moving = sitk.ReadImage("pet.nii.gz")
resampler = sitk.ResampleImageFilter()
resampler.SetReferenceImage(fixed) # 直接继承 Origin/Spacing/Direction/Size
resampler.SetInterpolator(sitk.sitkLinear) # PET 用线性,标签图用最近邻
resampler.SetDefaultPixelValue(0)
aligned = resampler.Execute(moving)
print("fixed origin:", fixed.GetOrigin())
print("aligned origin:", aligned.GetOrigin()) # 应与 fixed 一致
print("aligned direction:\n", aligned.GetDirection())
执行后应检查 aligned 的 Origin、Spacing、Direction 是否与 fixed 完全一致。若不一致,说明参考几何没有正确设置。对于标签图或分割掩膜,必须使用 sitkNearestNeighbor,否则插值会引入不存在的标签值。
几何一致性检查清单
在多图像处理流程中,建议在关键节点做一次几何核对。以下清单可直接用于工程检查:
| 检查项 | 方法 | 期望结果 |
|---|---|---|
| Origin 是否一致 | GetOrigin() 对比 | 对齐后与参考图相同 |
| Spacing 是否一致 | GetSpacing() 对比 | 对齐后与参考图相同 |
| Direction 是否一致 | GetDirection() 对比 | 对齐后与参考图相同 |
| Size 是否一致 | GetSize() 对比 | 对齐后与参考图相同 |
| 索引↔物理往返 | TransformIndexToPhysicalPoint 再转回 | 误差小于 1e-4 |
理解这四个几何量之后,后续的配准、重采样、ROI 提取、多模态融合才有可靠的空间基础。实际开发中,建议把几何检查封装成函数,在每次读取或生成图像后调用,避免因几何不一致导致的隐性错误。