2025-08-30 读取dicom CT Series，SimpleITK的使用

参考这里面的
医学图像处理 pydicom – 简书 (jianshu.com)

我觉得处理的时候，最好可以吧一个Series的CT slices单独放在一个文件夹下面。
读取一个系列的CT图形（CT Series），前面引文里的代码，先放着后面一段段看。

import os
import pydicom
import numpy
from matplotlib import pyplot

dcm = r D:code	emp3CT 
lstFilesDCM = []

# 将所有dicom文件读入
for diName, subdirList, fileList in os.walk(dcm):
    for filename in fileList:
        if ".dcm" in filename.lower():
            #print(filename)
            lstFilesDCM.append(os.path.join(diName, filename))

## 将第一张图片作为参考图
RefDs = pydicom.read_file(lstFilesDCM[10])
# print(RefDs)
# print(RefDs.pixel_array)
# print(RefDs.PatientPosition)
pyplot.imshow(RefDs.pixel_array, cmap=pyplot.cm.bone)
pyplot.show()

# 建立三维数组，长、宽、层数
ConstPixelDims = (int(RefDs.Rows), int(RefDs.Columns), len(lstFilesDCM))
print(ConstPixelDims)

# 得到 spacing 值 mm 单位
# PixelSpacing - 每个像素点实际的长度与宽度,单位(mm)
# SliceThickness - 每层切片的厚度,单位(mm)
ConstPixelSpacing = (float(RefDs.PixelSpacing[0]), float(RefDs.PixelSpacing[1]), float(RefDs.SliceThickness))

# 三维数据
# 0 到（第一个维数加一*像素间的间隔），步长为 ConstPixelSpacing
x = numpy.arange(0.0, (ConstPixelDims[0] + 1) * ConstPixelSpacing[0], ConstPixelSpacing[0])
y = numpy.arange(0.0, (ConstPixelDims[1] + 1) * ConstPixelSpacing[1], ConstPixelSpacing[1])
z = numpy.arange(0.0, (ConstPixelDims[2] + 1) * ConstPixelSpacing[2], ConstPixelSpacing[2])
print(len(x),"xxxx")

ArrayDicom = numpy.zeros(ConstPixelDims, dtype=RefDs.pixel_array.dtype)

# 遍历所有的 dicom 文件，读取图像数据，存放在 numpy 数组中
for filenameDCM in lstFilesDCM:
    ds = pydicom.read_file(filenameDCM)
    ArrayDicom[:, :, lstFilesDCM.index(filenameDCM)] = ds.pixel_array


# 轴状面显示
# 每英寸的像素数 dpi - dot per inch 越大图片越清晰。
pyplot.figure(dpi=1000)
# 将坐标轴都变为同等长度
# pyplot.axes().set_aspect( equal ,  datalim )
pyplot.axes().set_aspect( equal )
pyplot.set_cmap(pyplot.gray()) # 灰度化图片

# 第三个维度表明目前展示的是第几层 
pyplot.imshow(ArrayDicom[:, :, 66])
pyplot.show()

# 冠状面显示
pyplot.figure(dpi=100)
pyplot.axes().set_aspect( equal ,  datalim )
pyplot.set_cmap(pyplot.gray())
pyplot.imshow(ArrayDicom[:, 90, :])
pyplot.show()

lstFilesDCM是读取的CT slices文件名列表，貌似可以稍微改善一下:

确保所有读取的DICOM文件属于同一个序列并且只包含CT文件，您需要检查SeriesInstanceUID属性以及确保模态（Modality）是CT，另外还使用了stop_before_pixels=True参数来仅读取DICOM的元数据，这可以加速加载过程。如果读取过程中发生错误（例如，由于文件不是一个有效的DICOM文件），代码将简单地跳过该文件。

import os
import pydicom

def load_sorted_dicoms(folder_path):
    """
    Load and sort DICOMs in a folder based on SliceLocation and ensure they are all from the same CT series.
    """
    # List all files in the directory
    dicom_files = [os.path.join(folder_path, f) for f in os.listdir(folder_path) if os.path.isfile(os.path.join(folder_path, f))]

    dicoms = []
    for f in dicom_files:
        try:
            dcm = pydicom.dcmread(f, stop_before_pixels=True)  # Read the metadata only, to make it faster
            if hasattr(dcm,  Modality ) and dcm.Modality ==  CT :
                dicoms.append(dcm)
        except pydicom.errors.InvalidDicomError:
            # This file isn t a valid DICOM, just skip it
            pass

    # Ensure all DICOMs are from the same series
    series_uids = {dcm.SeriesInstanceUID for dcm in dicoms}
    if len(series_uids) != 1:
        raise ValueError("More than one DICOM series detected in the folder")

    # Sort the DICOMs based on the SliceLocation attribute
    dicoms.sort(key=lambda x: float(x.SliceLocation))

    return dicoms

folder_path =  path_to_your_folder 
sorted_dicoms = load_sorted_dicoms(folder_path)

或者可以使用SimpleITK来读取

import SimpleITK as sitk

def load_dicom_series(directory):
    # Get the directory names for all DICOM series in the given directory
    series_ids = sitk.ImageSeriesReader.GetGDCMSeriesIDs(directory)
    if not series_ids:
        raise ValueError("ERROR: given directory "" + directory + "" does not contain a DICOM series.")
    elif len(series_ids) > 1:
         raise ValueError(f"More than One series found in {directory}.")
    # For this example, we ll just take the first series. 
    # If there are multiple series, you ll need to handle them.
    series_file_names = sitk.ImageSeriesReader.GetGDCMSeriesFileNames(directory, series_ids[0])
    
    # Read the series
    series_reader = sitk.ImageSeriesReader()
    series_reader.SetFileNames(series_file_names)
    
    # Get the image
    image = series_reader.Execute()
    
    return image

directory = r D:code	emp3Test004615P_SMLC_45Gy_25F_(4956)CT 
dicom_image = load_dicom_series(directory)

精简版如下：

import SimpleITK as sitk
def load_dicom_series(directory):
    """加载DICOM序列到SimpleITK Image对象中"""
    reader = sitk.ImageSeriesReader()
    dicom_series = reader.GetGDCMSeriesFileNames(directory)
    reader.SetFileNames(dicom_series)
    image = reader.Execute()
    return image

获取一个slice作为参考图，读出三维方向的长宽高ConstPixelDims，建立一个numpy array把每个slice的pixel_array保存进去(貌似应该是排序好dcm.SliceLocation的)：

ConstPixelDims = (int(RefDs.Rows), int(RefDs.Columns), len(lstFilesDCM))
ArrayDicom = numpy.zeros(ConstPixelDims, dtype=RefDs.pixel_array.dtype)
for filenameDCM in lstFilesDCM:
    ds = pydicom.read_file(filenameDCM)
    ArrayDicom[:, :, lstFilesDCM.index(filenameDCM)] = ds.pixel_array
pyplot.imshow(ArrayDicom[:,:,60], cmap=pyplot.cm.bone)
pyplot.show()

使用 SimpleITK，我们不必手动排序DICOM文件或处理像素间距和切片厚度；它会为我们处理这些。
SimpleITK 读取DICOM序列后，我们可以很容易地将其转换为NumPy数组进行进一步的处理或可视化。

更新自己的评论：

用pyplot.imshow(ArrayDicom[:, :, 66])， 或者说66层的ds.pixel_array的话，正常显示，那就是说ds.pixel_array里面的维度也就是ndarray的(y, x)。 然后ds.Rows 行数 就是y方向的尺度，也就是Height图像高度，ds.Columns就是x方向的尺度，也就是图像的Width宽度

前面代码里的x、y、z

x, y, z数组是为了根据DICOM文件的元数据（即PixelSpacing和SliceThickness）定义三维空间的坐标系。具体来说：
x定义了图像的宽度方向上的坐标值。
y定义了图像的高度方向上的坐标值。
z定义了序列中切片的坐标值，即沿着深度方向（一般对应于病人从头到脚的方向）。
这些坐标系的定义有助于将像素空间转换为真实的物理空间（即毫米）。例如，如果一个DICOM图像的PixelSpacing属性是[0.5, 0.5]，这意味着每个像素对应于0.5毫米的真实空间。SliceThickness则是两个相邻切片之间的距离。