Python如何处理医疗数据？DICOM文件读取教程

pip install pydicom

import pydicom
import matplotlib.pyplot as plt
import numpy as np

try:
    # 读取DICOM文件
    ds = pydicom.dcmread("my_dicom_file.dcm")

    # 访问元数据，比如病人姓名
    patient_name = ds.PatientName
    modality = ds.Modality
    print(f"病人姓名: {patient_name}")
    print(f"影像模态: {modality}")

    # 获取像素数据
    # DICOM像素数据通常以原始整数形式存储，需要进行窗口化（Windowing）才能正确显示
    # 但为了简单起见，这里直接显示，实际应用中需要考虑WindowCenter和WindowWidth
    pixel_array = ds.pixel_array

    # 简单可视化，如果图像是灰度图
    if pixel_array.ndim == 2:
        plt.imshow(pixel_array, cmap=plt.cm.bone)
        plt.title(f"DICOM图像 - {patient_name} ({modality})")
        plt.axis('off')
        plt.show()
    elif pixel_array.ndim == 3: # 可能是彩色图或多帧图像
        # 对于彩色图 (Hounsfield units or RGB)
        # 需要根据PhotometricInterpretation进行处理
        # 这里仅作示例，直接显示第一帧或RGB图像
        if ds.PhotometricInterpretation == "RGB":
            plt.imshow(pixel_array)
            plt.title(f"DICOM图像 (RGB) - {patient_name} ({modality})")
            plt.axis('off')
            plt.show()
        else: # 可能是多帧灰度图
            print("这是多帧图像，显示第一帧。")
            plt.imshow(pixel_array[0], cmap=plt.cm.bone)
            plt.title(f"DICOM图像 (多帧) - {patient_name} ({modality})")
            plt.axis('off')
            plt.show()

except FileNotFoundError:
    print("错误：DICOM文件未找到，请确保路径正确。")
except Exception as e:
    print(f"处理DICOM文件时发生错误: {e}")

# 假设ds是已加载的DICOM数据集
try:
    print(f"病人ID: {ds.PatientID}")
    print(f"像素间距: {ds.PixelSpacing} mm")
    print(f"层厚: {ds.SliceThickness} mm")
    print(f"窗位/窗宽: {ds.WindowCenter} / {ds.WindowWidth}")
    print(f"图像位置: {ds.ImagePositionPatient}")
    print(f"图像方向: {ds.ImageOrientationPatient}")
except AttributeError as e:
    print(f"某个DICOM标签可能不存在: {e}")

import pydicom
import os
import numpy as np
import matplotlib.pyplot as plt

def load_dicom_series(dicom_dir):
    """
    加载一个目录下的所有DICOM文件，并尝试构建3D序列。
    """
    slices = []
    # 遍历目录，找到所有.dcm文件
    for s in os.listdir(dicom_dir):
        filepath = os.path.join(dicom_dir, s)
        if os.path.isfile(filepath) and s.endswith('.dcm'):
            try:
                ds = pydicom.dcmread(filepath)
                # 确保是图像文件，而不是其他DICOM对象（如结构报告等）
                if 'PixelData' in ds:
                    slices.append(ds)
            except Exception as e:
                print(f"跳过文件 {s}，因为它不是有效的DICOM图像文件或读取失败: {e}")

    # 根据InstanceNumber或ImagePositionPatient排序
    # InstanceNumber通常是简单的序列号，但ImagePositionPatient更可靠地反映物理位置
    # 这里我们优先使用ImagePositionPatient的Z轴分量进行排序
    if slices:
        # 尝试根据ImagePositionPatient的Z轴分量排序
        # 注意：ImagePositionPatient[2]是Z轴坐标
        slices.sort(key=lambda s: s.ImagePositionPatient[2] if 'ImagePositionPatient' in s else s.InstanceNumber)

        # 提取像素数据并堆叠
        # 注意原始像素值可能需要进行RescaleSlope和RescaleIntercept转换
        # 这里为了简化，直接获取pixel_array
        image_3d = np.stack([s.pixel_array for s in slices])

        # 获取体素间距，这对于3D可视化和测量很重要
        # 通常是 [SliceThickness, PixelSpacing[0], PixelSpacing[1]]
        try:
            slice_thickness = slices[0].SliceThickness if 'SliceThickness' in slices[0] else 1.0 # 默认1.0
            pixel_spacing = slices[0].PixelSpacing if 'PixelSpacing' in slices[0] else [1.0, 1.0] # 默认1.0
            voxel_spacing = [slice_thickness, pixel_spacing[0], pixel_spacing[1]]
        except AttributeError:
            print("警告: 无法获取完整的体素间距信息，使用默认值。")
            voxel_spacing = [1.0, 1.0, 1.0]

        return image_3d, voxel_spacing
    else:
        return None, None

# 示例用法
# 假设你的DICOM序列文件都在 'path/to/your/dicom_series/' 目录下
dicom_directory = 'path/to/your/dicom_series/' # 请替换为你的实际路径

if os.path.exists(dicom_directory):
    volume, spacing = load_dicom_series(dicom_directory)

    if volume is not None:
        print(f"成功加载3D体数据，形状: {volume.shape}")
        print(f"体素间距 (Z, Y, X): {spacing}")

        # 简单可视化：显示中间切片
        middle_slice_index = volume.shape[0] // 2
        plt.imshow(volume[middle_slice_index], cmap=plt.cm.bone)
        plt.title(f"3D序列中间切片 (Z={middle_slice_index})")
        plt.axis('off')
        plt.show()

        # 你现在可以对这个 `volume` (NumPy数组) 进行各种图像处理和分析了，
        # 比如阈值分割、边缘检测、3D渲染等。
    else:
        print("未找到有效的DICOM序列文件。")
else:
    print(f"错误: 目录 '{dicom_directory}' 不存在。")