RTSTRUCT文件3D结构Dice相似性指数(DSC)计算问题排查

遇到的问题 #

• 坐标提取索引越界：选择最后两个结构时，rtstruct.ROIContourSequence[roi_number -1]触发索引越界，强行调整索引为roi_number-10会导致结果错误。
• DSC结果无对称性：交换输入结构顺序，DSC从0.2853变为1.7147，违反DSC与输入顺序无关的特性。
• 修改坐标提取后DSC恒为1.0：调整坐标提取逻辑后，无论输入结构如何，DSC结果始终是1.0。
• Shapely多边形计算报错：提取的坐标数组仅含单个点，触发ValueError: A linearring requires at least 4 coordinates错误。

问题1：坐标提取索引越界的修复 #

问题根源 #

ROINumber是DICOM标签中的自定义编号，不一定和ROIContourSequence的列表索引一一对应（比如ROINumber可能从10开始编号，而列表索引从0开始），直接用roi_number-1作为索引必然会越界。

修复代码 #

def extract_coordinates(rtstruct, target_name):
    # 建立ROI编号到名称的映射
    roi_num_to_name = {roi.ROINumber: roi.ROIName.strip() for roi in rtstruct.StructureSetROISequence}
    # 遍历ROIContourSequence，通过ReferencedROINumber匹配目标结构
    for contour in rtstruct.ROIContourSequence:
        ref_roi_num = contour.ReferencedROINumber
        if roi_num_to_name.get(ref_roi_num, "").lower() == target_name.lower():
            coordinate_array = []
            for dataset in contour.ContourSequence:
                coords = dataset.ContourData
                # 按3个一组提取xyz坐标
                for i in range(0, len(coords), 3):
                    x, y, z = map(float, coords[i:i+3])
                    coordinate_array.append((x, y, z))
            return coordinate_array
    return None

问题2：DSC结果无对称性的修复 #

问题根源 #

初始代码错误使用np.logical_and和np.union1d处理Python集合，这两个函数仅适用于numpy数组，导致交集、并集计算逻辑错误，破坏了DSC的对称性。

修复代码 #

def calculate_dice_similarity(coords1, coords2):
    set1 = set(coords1)
    set2 = set(coords2)
    intersection = len(set1 & set2)
    # 严格遵循DSC公式，确保输入顺序不影响结果
    dice = (2.0 * intersection) / (len(set1) + len(set2))
    # 限制结果在0-1区间，避免浮点误差导致的微小超出
    return max(0.0, min(1.0, dice))

问题3：修改坐标提取后DSC恒为1.0的修复 #

问题根源 #

修改后的坐标提取函数中，return coordinate_array被放在了内层循环（for i in range(...)）内，导致仅提取第一个坐标点就直接返回，两个结构的坐标集合完全相同（仅单个点），因此DSC恒为1.0。

修复代码 #

将return coordinate_array移至所有循环外部，确保提取完整的坐标集合：

def extract_coordinates(rtstruct, target_name):
    roi_num_to_name = {roi.ROINumber: roi.ROIName.strip() for roi in rtstruct.StructureSetROISequence}
    for contour in rtstruct.ROIContourSequence:
        ref_roi_num = contour.ReferencedROINumber
        if roi_num_to_name.get(ref_roi_num, "").lower() == target_name.lower():
            coordinate_array = []
            for dataset in contour.ContourSequence:
                coords = dataset.ContourData
                for i in range(0, len(coords), 3):
                    x, y, z = map(float, coords[i:i+3])
                    coordinate_array.append((x, y, z))
            # 移至此处，遍历完所有轮廓点后再返回
            return coordinate_array
    return None

问题4：Shapely多边形计算报错的修复 #

问题根源 #

• Shapely的Polygon仅支持2D坐标，无法直接处理3D坐标。
• 单个点无法构成多边形，需按z轴切片分组，每个z层的轮廓点单独构成2D多边形。

修复代码 #

from shapely.geometry import Polygon
from collections import defaultdict

def group_coords_by_z(coords):
    # 按z坐标分组（做精度处理避免浮点误差）
    z_groups = defaultdict(list)
    for x, y, z in coords:
        z_key = round(z, 2)
        z_groups[z_key].append((x, y))
    return z_groups

def calculate_dice_shapely(coords1, coords2):
    z1 = group_coords_by_z(coords1)
    z2 = group_coords_by_z(coords2)
    all_z = set(z1.keys()).union(set(z2.keys()))
    
    total_intersection = 0
    total_area1 = 0
    total_area2 = 0
    
    for z in all_z:
        pts1 = z1.get(z, [])
        pts2 = z2.get(z, [])
        
        # 跳过无法构成有效多边形的情况
        if len(pts1) < 4 or len(pts2) < 4:
            continue
            
        try:
            poly1 = Polygon(pts1)
            poly2 = Polygon(pts2)
            
            area1 = poly1.area
            area2 = poly2.area
            intersection = poly1.intersection(poly2).area
            
            total_area1 += area1
            total_area2 += area2
            total_intersection += intersection
        except Exception:
            # 跳过自相交等无效多边形
            continue
    
    if total_area1 + total_area2 == 0:
        return 0.0
    # 按体积加权计算整体DSC
    return (2.0 * total_intersection) / (total_area1 + total_area2)

内容的提问来源于stack exchange，提问作者Francesco La Fauci