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 往往不再是单位阵。

判断一张图是否“轴对齐”,最直接的方法是检查 GetDirection() 是否接近单位矩阵。若明显偏离,说明索引轴与物理轴存在旋转关系,不能简单用索引差乘以 Spacing 来估算物理距离。

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 在映射中的位置。图中同时标出一个示例体素从索引到物理坐标的转换路径,便于对照理解。

索引空间 ↔ 物理空间的双向映射 索引空间 (i, j, k) 体素 (100, 80, 20) 离散整数索引,访问像素用 物理空间 (x, y, z) mm z x y p = (-123.4, -87.2, 45.6) 连续物理坐标,配准/标注用 仿射映射 Origin Spacing Direction 3×3 方向余弦 TransformIndexToPhysicalPoint TransformPhysicalPointToContinuousIndex p = Origin + Direction · (Spacing ⊙ (i, j, k)) | 索引 (0,0,0) 对应第一个体素中心

坐标系对齐的实操方法

当需要把一张图像对齐到另一张参考图像时,最稳妥的做法是显式构造目标几何,而不是手动改 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,否则插值会引入不存在的标签值。

常见错误:直接修改 GetOrigin() 返回的元组并不会改变图像几何,必须使用 SetOrigin()。另外,Origin 是索引 (0,0,0) 的物理坐标,不是图像包围盒的最小角点;当 Direction 含负方向时,两者可能相差很大。

几何一致性检查清单

在多图像处理流程中,建议在关键节点做一次几何核对。以下清单可直接用于工程检查:

检查项 方法 期望结果
Origin 是否一致 GetOrigin() 对比 对齐后与参考图相同
Spacing 是否一致 GetSpacing() 对比 对齐后与参考图相同
Direction 是否一致 GetDirection() 对比 对齐后与参考图相同
Size 是否一致 GetSize() 对比 对齐后与参考图相同
索引↔物理往返 TransformIndexToPhysicalPoint 再转回 误差小于 1e-4

理解这四个几何量之后,后续的配准、重采样、ROI 提取、多模态融合才有可靠的空间基础。实际开发中,建议把几何检查封装成函数,在每次读取或生成图像后调用,避免因几何不一致导致的隐性错误。