import pydicom
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.widgets import Slider
import os
from matplotlib.path import Path
import matplotlib.patches as patches
from matplotlib.colors import to_rgba
def load_ct_series(directory):
"""加载DICOM CT序列并提取切片位置信息"""
ct_slices = []
files = [f for f in os.listdir(directory) if f.endswith('.dcm')]
for file in files:
path = os.path.join(directory, file)
try:
ds = pydicom.dcmread(path)
if hasattr(ds, 'ImagePositionPatient'):
ct_slices.append(ds)
except Exception as e:
print(f"Error reading {file}: {e}")
# 按Z轴位置排序切片
ct_slices.sort(key=lambda ds: float(ds.ImagePositionPatient[2]))
return ct_slices
def load_rt_structure(rt_path, ct_series):
"""加载并解析RT Structure文件,匹配到CT切片"""
try:
rt_ds = pydicom.dcmread(rt_path)
except Exception as e:
print(f"Error reading RT Structure file: {e}")
return {}, {}
structures = {}
roi_info = {} # 存储所有ROI的信息,包括名称和颜色
# 提取所有ROI信息
rois = {}
if hasattr(rt_ds, 'StructureSetROISequence'):
for roi in rt_ds.StructureSetROISequence:
rois[roi.ROINumber] = roi.ROIName
else:
print("No StructureSetROISequence found in RT Structure file")
return {}, {}
# 提取轮廓序列
if hasattr(rt_ds, 'ROIContourSequence'):
for contour in rt_ds.ROIContourSequence:
roi_number = contour.ReferencedROINumber
if roi_number not in rois:
continue
roi_name = rois[roi_number]
color = contour.ROIDisplayColor if hasattr(contour, 'ROIDisplayColor') else [255, 0, 0]
# 存储ROI信息
roi_info[roi_name] = [c/255.0 for c in color] # 转换为0-1范围
# 提取轮廓数据
if hasattr(contour, 'ContourSequence'):
for contour_sequence in contour.ContourSequence:
# 提取轮廓点数据 (x,y,z)
contour_data = np.array(contour_sequence.ContourData)
points = contour_data.reshape((-1, 3))
# 获取轮廓的Z坐标(使用第一个点的Z坐标作为参考)
z_pos = points[0, 2]
# 找到最接近的CT切片
slice_index = find_closest_slice(ct_series, z_pos, tolerance=1.0) # 增加容差
if slice_index is not None:
if slice_index not in structures:
structures[slice_index] = []
# 存储轮廓信息
structures[slice_index].append({
'points': points[:, :2], # 只取x,y坐标
'roi_name': roi_name,
'color': [c/255.0 for c in color] # 转换为0-1范围
})
return structures, roi_info # 返回轮廓结构和ROI信息
def find_closest_slice(ct_series, z_pos, tolerance=1.0):
"""查找最接近给定Z位置的CT切片索引"""
min_distance = float('inf')
closest_index = None
for i, ds in enumerate(ct_series):
slice_z = float(ds.ImagePositionPatient[2])
distance = abs(slice_z - z_pos)
if distance < min_distance:
min_distance = distance
closest_index = i
# 检查是否在容差范围内
if min_distance <= tolerance:
return closest_index
else:
return None
def apply_window_level(image, window_center, window_width):
"""应用窗宽窗位调整"""
min_val = window_center - window_width / 2
max_val = window_center + window_width / 2
windowed = np.clip(image, min_val, max_val)
return (windowed - min_val) / (max_val - min_val)
def patient_coord_to_pixel_coord(patient_coord, ds):
"""将患者坐标系中的点转换到像素坐标系"""
# 获取DICOM图像的方向和位置信息
image_orientation = ds.ImageOrientationPatient
image_position = ds.ImagePositionPatient
pixel_spacing = ds.PixelSpacing
# 创建转换矩阵
x_vec = np.array(image_orientation[0:3])
y_vec = np.array(image_orientation[3:6])
z_vec = np.cross(x_vec, y_vec)
# 计算转换矩阵
transform_matrix = np.vstack((x_vec, y_vec, z_vec)).T
transform_matrix = transform_matrix * np.array([pixel_spacing[0], pixel_spacing[1], 1])
# 计算患者坐标到图像坐标的偏移
offset = np.array(patient_coord) - np.array(image_position)
# 应用转换
pixel_coord = np.linalg.solve(transform_matrix.T, offset.T)
return pixel_coord[0:2] # 返回x,y像素坐标
def update_display(slice_idx):
"""更新显示的CT切片和轮廓"""
global current_slice_idx, axial_img, contour_patches, text_annotations
# 清除之前的图像和轮廓
if axial_img is not None:
axial_img.remove()
for patch in contour_patches:
patch.remove()
for text in text_annotations:
text.remove()
contour_patches = []
text_annotations = []
# 获取当前切片
ds = ct_series[slice_idx]
pixel_array = ds.pixel_array
# 转换为HU单位(如果适用)
if hasattr(ds, 'RescaleSlope') and hasattr(ds, 'RescaleIntercept'):
intercept = ds.RescaleIntercept
slope = ds.RescaleSlope
pixel_array = pixel_array * slope + intercept
# 应用窗宽窗位(默认为软组织窗)
window_center = 40 if not hasattr(ds, 'WindowCenter') else ds.WindowCenter
window_width = 400 if not hasattr(ds, 'WindowWidth') else ds.WindowWidth
# 确保窗宽窗位是数值而不是序列
if isinstance(window_center, pydicom.multival.MultiValue):
window_center = window_center[0]
if isinstance(window_width, pydicom.multival.MultiValue):
window_width = window_width[0]
windowed_image = apply_window_level(pixel_array, window_center, window_width)
# 显示CT图像
axial_img = ax.imshow(windowed_image, cmap='gray')
ax.set_title(f'Slice {slice_idx+1}/{len(ct_series)}\nZ:{float(ds.ImagePositionPatient[2]):.2f} mm')
ax.axis('off')
# 如果有该切片的轮廓,则添加绘图
if slice_idx in rt_structures:
for contour in rt_structures[slice_idx]:
points = contour['points']
color = contour['color']
roi_name = contour['roi_name']
# 转换坐标到像素空间
pixel_points = []
for point in points:
# 添加Z坐标(使用当前切片的Z位置)
patient_coord = [point[0], point[1], float(ds.ImagePositionPatient[2])]
pixel_coord = patient_coord_to_pixel_coord(patient_coord, ds)
pixel_points.append(pixel_coord)
pixel_points = np.array(pixel_points)
# 转换为图像坐标
path = Path(pixel_points)
patch = patches.PathPatch(path, facecolor='none',
edgecolor=color, linewidth=2)
ax.add_patch(patch)
contour_patches.append(patch)
# 添加轮廓名称文本
if len(pixel_points) > 0:
# 计算轮廓中心点
center = np.mean(pixel_points, axis=0)
text = ax.text(center[0], center[1], roi_name,
color=color, fontsize=8, ha='center', va='center',
bbox=dict(boxstyle="round,pad=0.3", facecolor="black", alpha=0.5))
text_annotations.append(text)
plt.draw()
def on_slider_change(val):
"""滑块值改变时的回调函数"""
slice_idx = int(slice_slider.val)
update_display(slice_idx)
# 主程序
if __name__ == "__main__":
# 配置路径
ct_directory = '/Users/yangxiong/Desktop/lattice-RT/933022' # CT图像目录
rt_path = '/Users/yangxiong/Desktop/lattice-RT/933022/RS.1.3.46.670589.33.1.63893199939983227300001.4937323749018908579.1.dcm' # RT Structure文件路径
# 加载CT序列
print("Loading CT series...")
ct_series = load_ct_series(ct_directory)
print(f"Loaded {len(ct_series)} CT slices")
# 加载并处理RT Structure
print("Loading RT Structure...")
rt_structures, roi_info = load_rt_structure(rt_path, ct_series)
print(f"Found contours for {len(rt_structures)} slices")
print(f"Found {len(roi_info)} unique ROIs: {list(roi_info.keys())}")
# 创建图形和轴
fig, ax = plt.subplots(1, 1, figsize=(10, 10))
plt.subplots_adjust(bottom=0.2) # 为滑块留出空间
# 创建滑块
slider_ax = plt.axes([0.2, 0.05, 0.65, 0.03])
slice_slider = Slider(slider_ax, 'Slice', 0, len(ct_series)-1, valinit=0, valfmt='%0.0f')
# 注册滑块回调函数
slice_slider.on_changed(on_slider_change)
# 全局变量
current_slice_idx = 0
axial_img = None
contour_patches = []
text_annotations = []
# 创建图例
legend_handles =
legend_labels = []
# 为每个ROI创建一个图例句柄
for roi_name, color in roi_info.items():
patch = patches.Patch(color=color, label=roi_name)
legend_handles.append(patch)
legend_labels.append(roi_name)
# 添加图例
if legend_handles:
fig.legend(legend_handles, legend_labels, loc='upper right', bbox_to_anchor=(0.95, 0.95))
# 初始显示
update_display(0)
plt.show()修改代码使滑片能滑动