# -*- coding: utf-8 -*-
import os
import pydicom
import numpy as np
import matplotlib.pyplot as plt
from pydicom.errors import InvalidDicomError
def load_dicom_series(series_path):
"""
加载指定路径下的所有DICOM文件，并按照切片位置排序
参数:
series_path (str): DICOM系列所在的目录路径
返回:
list: 排序后的DICOM数据集列表
"""
# 1. 检查路径是否存在
if not os.path.exists(series_path):
raise FileNotFoundError(f"目录不存在: {series_path}")
# 2. 获取目录下所有DICOM文件
dicom_files = []
for root, _, files in os.walk(series_path):
for file in files:
file_path = os.path.join(root, file)
try:
# 3. 尝试读取DICOM文件
dataset = pydicom.dcmread(file_path)
# 4. 验证是否包含必要属性
if hasattr(dataset, 'SliceLocation'):
dicom_files.append(dataset)
else:
print(f"警告: 文件缺少SliceLocation属性，跳过: {file_path}")
except InvalidDicomError:
# 5. 跳过非DICOM文件
print(f"跳过非DICOM文件: {file_path}")
# 6. 检查是否找到DICOM文件
if not dicom_files:
raise ValueError("目录中未找到有效的DICOM文件")
# 7. 按切片位置排序
dicom_files.sort(key=lambda ds: float(ds.SliceLocation))
return dicom_files
def extract_pixel_data(dicom_series):
"""
从DICOM系列中提取像素数据并转换为3D数组
参数:
dicom_series (list): 排序后的DICOM数据集列表
返回:
tuple: (3D像素数组, 元数据字典)
"""
# 1. 创建3D数组存储所有切片
first_slice = dicom_series[0]
rows = first_slice.Rows
columns = first_slice.Columns
num_slices = len(dicom_series)
# 2. 初始化3D数组
volume = np.zeros((num_slices, rows, columns), dtype=np.int16)
# 3. 提取元数据
metadata = {
'patient_id': getattr(first_slice, 'PatientID', '未知'),
'patient_name': getattr(first_slice, 'PatientName', '未知'),
'study_date': getattr(first_slice, 'StudyDate', '未知'),
'modality': getattr(first_slice, 'Modality', '未知'),
'series_description': getattr(first_slice, 'SeriesDescription', '未知'),
'pixel_spacing': getattr(first_slice, 'PixelSpacing', [1.0, 1.0]),
'slice_thickness': getattr(first_slice, 'SliceThickness', 1.0),
'image_orientation': getattr(first_slice, 'ImageOrientationPatient', None),
'image_position': getattr(first_slice, 'ImagePositionPatient', None),
'window_center': getattr(first_slice, 'WindowCenter', None),
'window_width': getattr(first_slice, 'WindowWidth', None)
}
# 4. 处理每个切片
for i, ds in enumerate(dicom_series):
# 获取像素数组
pixel_array = ds.pixel_array
# 应用Rescale Slope和Intercept（如果存在）
if hasattr(ds, 'RescaleSlope') and hasattr(ds, 'RescaleIntercept'):
pixel_array = pixel_array * ds.RescaleSlope + ds.RescaleIntercept
# 存储到3D数组中
volume[i, :, :] = pixel_array
return volume, metadata
def visualize_dicom_series(volume, metadata, slice_index=None):
"""
可视化DICOM系列数据
参数:
volume (np.array): 3D像素数组
metadata (dict): DICOM元数据
slice_index (int): 指定要显示的切片索引（如为None则显示中间切片）
"""
# 确定要显示的切片
if slice_index is None:
slice_index = volume.shape[0] // 2
# 获取切片数据
slice_data = volume[slice_index, :, :]
# 创建图像
plt.figure(figsize=(10, 10))
# 设置窗口参数（如果存在）
if metadata['window_center'] and metadata['window_width']:
if isinstance(metadata['window_center'], pydicom.multival.MultiValue):
wc = float(metadata['window_center'][0])
ww = float(metadata['window_width'][0])
else:
wc = float(metadata['window_center'])
ww = float(metadata['window_width'])
vmin = wc - ww/2
vmax = wc + ww/2
plt.imshow(slice_data, cmap='gray', vmin=vmin, vmax=vmax)
else:
# 自动调整对比度
plt.imshow(slice_data, cmap='gray')
# 添加标题和信息
title = f"{metadata['modality']} 扫描\n"
title += f"患者: {metadata['patient_name']} (ID: {metadata['patient_id']})\n"
title += f"日期: {metadata['study_date']}, 系列: {metadata['series_description']}\n"
title += f"切片: {slice_index+1}/{volume.shape[0]} (位置: {slice_index * metadata['slice_thickness']:.2f} mm)"
plt.title(title, fontsize=12)
plt.axis('off')
plt.colorbar(label='Hounsfield Units (HU)' if metadata['modality'] == 'CT' else 'Pixel Value')
plt.show()
def main():
# 设置DICOM系列路径
series_path = "series-00001" # 替换为你的实际路径
try:
# 1. 加载DICOM系列
print(f"正在加载DICOM系列: {series_path}")
dicom_series = load_dicom_series(series_path)
print(f"成功加载 {len(dicom_series)} 个DICOM文件")
# 2. 提取像素数据和元数据
volume, metadata = extract_pixel_data(dicom_series)
print("\n元数据信息:")
print(f"患者ID: {metadata['patient_id']}")
print(f"患者姓名: {metadata['patient_name']}")
print(f"检查日期: {metadata['study_date']}")
print(f"模态: {metadata['modality']}")
print(f"系列描述: {metadata['series_description']}")
print(f"体素尺寸: {metadata['pixel_spacing'][0]}×{metadata['pixel_spacing'][1]} mm²")
print(f"切片厚度: {metadata['slice_thickness']} mm")
print(f"数据形状: {volume.shape} (切片×行×列)")
# 3. 可视化DICOM数据
print("\n可视化中间切片...")
visualize_dicom_series(volume, metadata)
# 4. 可选：保存为NIFTI格式（需要nibabel库）
# save_as_nifti(volume, metadata)
except Exception as e:
print(f"错误: {str(e)}")
if __name__ == "__main__":
main()
读取、处理和展示 DICOM 格式的医学影像文件：
1. 导入所需库
- os：用于文件和目录操作。
- pydicom：专门用于读取 DICOM 格式文件的库。
- numpy：用于数值计算和数组操作。
- matplotlib.pyplot：用于绘制图像和图表。
- InvalidDicomError：pydicom 库中的异常类，用于捕获无效的 DICOM 文件。
2. 函数 load_dicom_series
- 作用：从指定的文件夹中读取所有的 DICOM 文件，并按照切片位置进行排序。
- 步骤：
- 检查指定的文件夹是否存在，如果不存在则报错。
- 遍历文件夹及其子文件夹，读取每个文件。对于每个文件，尝试使用 pydicom.dcmread 读取 DICOM 文件。
- 检查读取的 DICOM 文件是否包含 SliceLocation 属性，确保文件是有效的医学影像切片文件。如果缺少该属性，会发出警告并跳过该文件。
- 如果读取文件时出现 InvalidDicomError 异常，说明该文件不是有效的 DICOM 文件，会跳过该文件。
- 最后，按照 SliceLocation 属性对读取的 DICOM 文件进行排序，并返回排序后的 DICOM 数据集列表。
3. 函数 extract_pixel_data
- 作用：从排序后的 DICOM 数据集列表中提取像素数据，并将其转换为一个 3D 数组。同时，提取相关的元数据信息。
- 步骤：
- 获取第一个切片的行数、列数和切片数量，初始化一个 3D 数组 volume 用于存储所有切片的像素数据。
- 从第一个切片中提取元数据信息，包括患者 ID、姓名、检查日期、模态（Modality）、系列描述、像素间距、切片厚度、图像方向、图像位置、窗宽和窗位等。
- 遍历每个切片，获取其像素数组，并应用 RescaleSlope 和 RescaleIntercept（如果存在）进行像素值的调整。最后将调整后的像素数组存储到 3D 数组 volume 中。
4. 函数 visualize_dicom_series
- 作用：可视化 DICOM 系列数据。
- 步骤：
- 如果未指定 slice_index，则显示中间的切片。
- 根据元数据中的窗宽（window_width）和窗位（window_center）设置图像的显示范围，以便更好地显示医学影像的细节。如果没有窗宽和窗位信息，则自动调整对比度。
- 添加标题和信息，包括模态、患者信息、检查日期、系列描述、当前显示的切片索引和位置等。
- 最后显示图像，并添加颜色条，用于表示像素值的范围。
5. 主函数 main
- 作用：程序的入口点。
- 步骤：
- 设置 DICOM 系列的路径 series_path，你需要将其替换为实际的 DICOM 文件所在的目录路径。
- 在 try 块中，依次调用 load_dicom_series 加载 DICOM 系列，调用 extract_pixel_data 提取像素数据和元数据，并调用 visualize_dicom_series 可视化 DICOM 数据。
- 如果在执行过程中出现异常，会捕获异常并打印错误信息。
这段代码提供了一个完整的流程，用于读取、处理和可视化 DICOM 格式的医学影像数据，方便医学影像的分析和诊断。
3D数组
1. 三个维度
第一个维度：
医学影像（如 CT、MRI 等）通常是由多个二维切片（slice）按一定顺序堆叠而成的三维结构
第二个维度（行数，rows）
- 含义：表示每个切片的高度，即切片图像在垂直方向上的像素数量。
- 业务需求：在医学影像中，每个切片都是一个二维图像，行数决定了图像的高度。例如，一个 CT 切片可能有 512 行像素，这意味着图像的高度为 512 个像素。
第三个维度（列数，columns）
- 含义：表示每个切片的宽度，即切片图像在水平方向上的像素数量。
- 业务需求：列数决定了图像的宽度。例如，一个 CT 切片可能有 512 列像素，这意味着图像的宽度为 512 个像素。
2. 空间关系的保持
每个切片在三维空间中都有其特定的位置（由 SliceLocation 属性表示）。使用 3D 数组可以保持这些切片之间的空间关系，使得在处理和分析时能够正确地理解和利用这种空间信息。例如，在进行三维重建或体积计算时，需要知道每个切片在三维空间中的相对位置。
3. 方便后续处理
许多医学影像处理任务（如三维重建、体积渲染、多平面重建等）都需要访问和处理整个三维数据。使用 3D 数组可以方便地进行这些操作。例如，通过对 3D 数组进行切片、旋转、缩放等操作，可以生成不同视角的图像或进行其他高级处理。
4. 内存效率
虽然 3D 数组可能占用较多的内存，但在大多数情况下，这种内存开销是可以接受的。而且，使用 3D 数组可以避免在处理过程中频繁地读取和写入磁盘，从而提高处理效率。
5. 与现有库和工具的兼容性
许多医学影像处理库和工具（如 ITK、VTK、SimpleITK 等）都支持 3D 数组作为输入和输出格式。使用 3D 数组可以方便地与这些库和工具进行集成，从而利用它们提供的高级功能。
扫一扫
1157
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
