3. 图像基础操作与像素访问

在 SimpleITK 中，图像（sitk.Image）不仅是一个多维数组，更是一个带有物理空间信息的对象。它同时维护两套坐标体系：索引坐标（Index）以整数描述像素在数组中的位置，物理坐标（Physical Point）以毫米等物理单位描述该像素在真实空间中的位置。理解这两套坐标的转换关系，是后续配准、重采样、ROI 提取等操作的基础。本节围绕像素读写、坐标转换、信息查询与统计四类操作展开。

一、像素访问：GetPixel 与 SetPixel

SimpleITK 提供 GetPixel 与 SetPixel 两个方法，分别用于读取和写入单个像素值。它们接受索引坐标，返回或设置的是像素的标量值（对于标量图像）或像素 ID（对于标签图像）。

import SimpleITK as sitk

img = sitk.ReadImage("ct_volume.nii.gz")

# 读取索引 (x=100, y=80, z=20) 处的像素值
value = img.GetPixel(100, 80, 20)
print("像素值:", value)

# 写入像素值
img.SetPixel(100, 80, 20, 1024)

需要注意：GetPixel 的索引顺序是 (x, y, z)，与 NumPy 数组的 [z, y, x] 顺序相反。这是初学者最容易出错的地方。若图像是向量图像（如 RGB 或位移场），GetPixel 返回的是长度为像素分量数的元组。

对于批量访问，逐像素调用 GetPixel 效率很低。此时应使用 GetArrayFromImage() 将图像转为 NumPy 数组，在数组上做向量化操作，再用 GetImageFromArray() 转回。但要注意：转回时若未显式设置原点、间距和方向，物理空间信息会丢失。

二、索引与物理坐标转换

SimpleITK 提供四个核心方法完成索引与物理坐标的互转：

• TransformIndexToPhysicalPoint(index)：索引 → 物理坐标
• TransformPhysicalPointToIndex(point)：物理坐标 → 索引
• TransformPhysicalPointToContinuousIndex(point)：物理坐标 → 连续索引（浮点）
• TransformContinuousIndexToPhysicalPoint(cindex)：连续索引 → 物理坐标

物理坐标由原点（Origin）、间距（Spacing）和方向（Direction）共同决定。对于二维图像，转换公式为：

P = O + D · (S ⊙ I)

其中：
  P = 物理坐标 (px, py)
  O = 原点 (ox, oy)
  D = 方向矩阵（2×2）
  S = 间距 (sx, sy)
  I = 索引 (ix, iy)
  ⊙ 表示逐元素相乘

下面是一个具体示例。假设某 CT 图像：

origin  = (10.0, 20.0, 30.0)   # mm
spacing = (0.7, 0.7, 1.5)      # mm
direction = 单位矩阵

index = (100, 80, 20)
physical = img.TransformIndexToPhysicalPoint(index)
# 计算：x = 10 + 0.7*100 = 80.0
#       y = 20 + 0.7*80  = 76.0
#       z = 30 + 1.5*20  = 60.0
# 结果：(80.0, 76.0, 60.0)

反向转换时，若物理坐标不落在整数索引上，TransformPhysicalPointToIndex 会返回最近的整数索引，而 TransformPhysicalPointToContinuousIndex 返回浮点索引，便于插值。在配准和重采样中，通常使用连续索引配合插值器获取亚像素精度。

三、图像信息查询

SimpleITK 提供一组查询方法，用于获取图像的元数据：

这些信息在读写图像、重采样、配准时必须保持一致。例如，将 NumPy 数组转回 SimpleITK 图像时，应显式调用 CopyInformation() 或手动设置上述元数据，否则物理坐标信息会丢失。

四、图像统计

SimpleITK 提供 StatisticsImageFilter 用于计算图像的全局统计量，包括最小值、最大值、均值、方差、总和等。它支持在指定区域内计算，也支持按标签统计。

stats = sitk.StatisticsImageFilter()
stats.Execute(img)

print("最小值:", stats.GetMinimum())
print("最大值:", stats.GetMaximum())
print("均值:", stats.GetMean())
print("方差:", stats.GetVariance())
print("总和:", stats.GetSum())

若需按标签区域统计，可使用 LabelStatisticsImageFilter，它同时接受灰度图像和标签图像，输出每个标签对应的统计量。这在肿瘤体积测量、器官分割评估中非常常用。

label_stats = sitk.LabelStatisticsImageFilter()
label_stats.Execute(img, label_img)

for label in label_stats.GetLabels():
    print(f"标签 {label}: 均值={label_stats.GetMean(label):.2f}, "
          f"体素数={label_stats.GetCount(label)}")

统计结果可用于判断图像是否正常。例如，CT 图像中空气的 HU 值约为 -1000，若最小值远高于此，可能说明图像被截断或窗宽窗位设置异常。均值与方差则可用于评估图像对比度和噪声水平。

五、综合示例：ROI 提取与统计

下面给出一个完整示例：在物理坐标 (80.0, 76.0, 60.0) 处提取 32×32×8 的 ROI，并计算其统计量。

import SimpleITK as sitk
import numpy as np

img = sitk.ReadImage("ct_volume.nii.gz")

# 物理坐标转索引
center = (80.0, 76.0, 60.0)
idx = img.TransformPhysicalPointToIndex(center)

# 定义 ROI 大小和起始索引
size = (32, 32, 8)
start = (idx[0] - 16, idx[1] - 16, idx[2] - 4)

# 提取 ROI
roi = sitk.RegionOfInterest(img, size, start)

# 统计
stats = sitk.StatisticsImageFilter()
stats.Execute(roi)
print("ROI 均值:", stats.GetMean())
print("ROI 标准差:", np.sqrt(stats.GetVariance()))

该示例展示了索引与物理坐标转换、ROI 提取和统计的联合使用。实际工程中，ROI 的起始索引需做边界检查，避免越界。SimpleITK 在越界时会抛出异常，因此建议在提取前用 GetSize() 验证范围。

掌握上述四类操作后，即可在 SimpleITK 中完成从像素级读写到区域统计的完整流程。实际开发中，建议始终先查询图像元数据，确认索引与物理坐标的对应关系，再进行后续处理，避免因坐标顺序或方向矩阵理解错误导致结果偏差。