SimpleITK 获取CT spacing 底层原理

SimpleITK 获取CT spacing 底层原理

一、层厚、层间距概念

层厚： CT扫描机扫描出来的断层的层的厚度， 通常用Slice thickness表示， 比如 5mm

层间隔：一般用 Slice interval 或 Slice increment来表示，比如 5mm。在 CT 中，层间隔表示相邻两层中心之间的距离

如上图所示， 在CT中， 相邻两层CT中心之间的距离代表层间距， 当层厚==层间距时，无重叠， 代表CT断层之间是连续的， 那么同理， 当层厚！=层间距时， 代表断层之间是有距离或者重叠的

二、pydicom中的层厚以及spacing

spacing:

一张二维的医学图像，图像尺寸为[6x7]，每个像素之间的物理空间距离（spacing），x轴、y轴分别是（20，30）。扩展到三维医学图像，就出现【体素】，就是三维的“像素”。相比像素有三个轴（x、y、z），同样每个轴都有spacing。 意思相同大小尺寸的医学图片，spacing也是不同。

如下是pydicom以及SimpleITK获取spacing：

import pydicom as dicom
import numpy as np
from pathlib import Path

PathDicom = 'dcm序列所在文件夹路径'

data_dir = Path(PathDicom)
dicom_files = list(data_dir.glob('*.dcm')) # 获取目录下所有 .dcm 文件

dicom_files.sort(key=lambda x: int(x.stem.rsplit('.', 1)[1])) #将 .dcm文件 按序号排列 注意： sorted() 方法会按字符排列， 此处需要按末尾数字排列
slices = [dicom.read_file(dicom_file, force=True) for dicom_file in dicom_files]

# 采用以下方法返回的SliceThickness是错误的
print(slices[0].SliceThickness)

# 正确的层厚获取方法
thickness = np.abs(slices[0].SliceLocation - slices[1].SliceLocation) # 获取 dcm 文件的层厚

# 获取整个 CT 序列的 spacing, 等同于 SimpleITK 的 GetSpacing() 方法
spacing_z_list = []
for i in range(1, len(slices)):
    spacing_z_list.append(slices[i].SliceLocation - slices[i - 1].SliceLocation)


spacing_z = abs(sum(spacing_z_list) / len(spacing_z_list)) # 获取Z轴的spacing
print(f"spacing_z : {spacing_z}")

如上代码所示， 因为dicom文件层次之间可能会存在重叠或者间隔， 故z轴的spacing应该取所有相邻序列层厚的均值， 即abs(sum(spacing_z_list) / len(spacing_z_list))。

注意：在临床CT数据中， 有时候会遇到一些数据在ITK-SNAP软件中显示如下， 有多个component：

这种情况下， 使用SimpleITK获取的Z轴spacing就是利用我上述代码的逻辑计算的， 并且单独的component之间数据的层厚是相同的， 只是因为两个component连接处的层厚大大超出单独component内的层厚， 故不注意的话会造成一些问题。这个component信息可以通过pydicom包中AcquisitionNumber获取

singleDicomFile # 单独的一个dcm文件
print(singleDicomFile.AcquisitionNumber) # 查看某个dcm文件所属component

参考链接：

【CTMR技术】CT与MR中的层厚&层间隔一样吗？

DICOM数据字段说明

关于DICOM文件的SliceThickness参数错误问题

【碎碎念】医学图像处理工具SimpleITK与DataSet