欢迎您访问程序员文章站本站旨在为大家提供分享程序员计算机编程知识!
您现在的位置是: 首页

【SimpleITK】医疗影像分割结果评价指标计算

程序员文章站 2022-05-08 14:37:53
...

Overlap Measures

  • jaccard
  • dice
  • volume_similarity
  • false_negative
  • false_positive

其中:
Volume measures:
volumeSimilarity2(v1v2)v1+v2 volumeSimilarity \frac{2*(v1-v2)}{v1+v2}

建立枚举对象:

from enum import Enum
# Use enumerations to represent the various evaluation measures
class OverlapMeasures(Enum):
    jaccard, dice, volume_similarity, false_negative, false_positive = range(5)  
reference_segmentation = reference_segmentation_STAPLE

建立空的数组,为了后面保存结果:

# Empty numpy arrays to hold the results 
overlap_results = np.zeros((len(segmentations),len(OverlapMeasures.__members__.items())))  
overlap_results

array([[0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0.]])

计算每个指标:

overlap_measures_filter = sitk.LabelOverlapMeasuresImageFilter()
for i, seg in enumerate(segmentations):
    # Overlap measures
    overlap_measures_filter.Execute(reference_segmentation, seg)
    overlap_results[i,OverlapMeasures.jaccard.value] = overlap_measures_filter.GetJaccardCoefficient()
    overlap_results[i,OverlapMeasures.dice.value] = overlap_measures_filter.GetDiceCoefficient()
    overlap_results[i,OverlapMeasures.volume_similarity.value] = overlap_measures_filter.GetVolumeSimilarity()
    overlap_results[i,OverlapMeasures.false_negative.value] = overlap_measures_filter.GetFalseNegativeError()
    overlap_results[i,OverlapMeasures.false_positive.value] = overlap_measures_filter.GetFalsePositiveError()
overlap_results

结果:

array([[ 0.82 , 0.901, 0.052, 0.075, 0.122],
[ 0.881, 0.937, -0.013, 0.069, 0.057],
[ 0.843, 0.915, -0.088, 0.124, 0.044]])

Surface Distance Measures

  • hausdorff_distance
  • mean_surface_distance
  • median_surface_distance
  • std_surface_distance
  • max_surface_distance

建立枚举对象:

class SurfaceDistanceMeasures(Enum):
    hausdorff_distance, mean_surface_distance, median_surface_distance, std_surface_distance, max_surface_distance = range(5)   
reference_segmentation = reference_segmentation_STAPLE

建立空的数组,为了后面保存结果:

surface_distance_results = np.zeros((len(segmentations),len(SurfaceDistanceMeasures.__members__.items()))) 
surface_distance_results

array([[0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0.],
[0., 0., 0., 0., 0.]])

对GT进行预处理:

# Use the absolute values of the distance map to compute the surface distances (distance map sign, outside or inside 
# relationship, is irrelevant)
label = 1
reference_distance_map = sitk.Abs(sitk.SignedMaurerDistanceMap(reference_segmentation, squaredDistance=False))
reference_surface = sitk.LabelContour(reference_segmentation)

statistics_image_filter = sitk.StatisticsImageFilter()
# Get the number of pixels in the reference surface by counting all pixels that are 1.
statistics_image_filter.Execute(reference_surface)
num_reference_surface_pixels = int(statistics_image_filter.GetSum()) 

计算指标:

hausdorff_distance_filter = sitk.HausdorffDistanceImageFilter()
for i, seg in enumerate(segmentations):
    hausdorff_distance_filter.Execute(reference_segmentation, seg)
    surface_distance_results[i,SurfaceDistanceMeasures.hausdorff_distance.value] = hausdorff_distance_filter.GetHausdorffDistance()
    # Symmetric surface distance measures
    segmented_distance_map = sitk.Abs(sitk.SignedMaurerDistanceMap(seg, squaredDistance=False, useImageSpacing=True))
    segmented_surface = sitk.LabelContour(seg)
        
    # Multiply the binary surface segmentations with the distance maps. The resulting distance
    # maps contain non-zero values only on the surface (they can also contain zero on the surface)
    seg2ref_distance_map = reference_distance_map*sitk.Cast(segmented_surface, sitk.sitkFloat32)
    ref2seg_distance_map = segmented_distance_map*sitk.Cast(reference_surface, sitk.sitkFloat32)
        
    # Get the number of pixels in the reference surface by counting all pixels that are 1.
    statistics_image_filter.Execute(segmented_surface)
    num_segmented_surface_pixels = int(statistics_image_filter.GetSum())
    
    # Get all non-zero distances and then add zero distances if required.
    seg2ref_distance_map_arr = sitk.GetArrayViewFromImage(seg2ref_distance_map)
    seg2ref_distances = list(seg2ref_distance_map_arr[seg2ref_distance_map_arr!=0]) 
    seg2ref_distances = seg2ref_distances + \
                        list(np.zeros(num_segmented_surface_pixels - len(seg2ref_distances)))
    ref2seg_distance_map_arr = sitk.GetArrayViewFromImage(ref2seg_distance_map)
    ref2seg_distances = list(ref2seg_distance_map_arr[ref2seg_distance_map_arr!=0]) 
    ref2seg_distances = ref2seg_distances + \
                        list(np.zeros(num_reference_surface_pixels - len(ref2seg_distances)))
        
    all_surface_distances = seg2ref_distances + ref2seg_distances

    # The maximum of the symmetric surface distances is the Hausdorff distance between the surfaces. In 
    # general, it is not equal to the Hausdorff distance between all voxel/pixel points of the two 
    # segmentations, though in our case it is. More on this below.
    surface_distance_results[i,SurfaceDistanceMeasures.mean_surface_distance.value] = np.mean(all_surface_distances)
    surface_distance_results[i,SurfaceDistanceMeasures.median_surface_distance.value] = np.median(all_surface_distances)
    surface_distance_results[i,SurfaceDistanceMeasures.std_surface_distance.value] = np.std(all_surface_distances)
    surface_distance_results[i,SurfaceDistanceMeasures.max_surface_distance.value] = np.max(all_surface_distances)
print(surface_distance_results)

结果:

[[4.905 0.459 0. 0.85 4.905]
[3.469 0.29 0. 0.689 3.469]
[5.203 0.431 0. 0.831 5.203]]

some API

  1. sitk.LabelContour
    Labels the pixels on the border of the objects in a labeled image. 得到mask的3维轮廓。
    对reference_segmentation进行轮廓提取处理:
reference_surface = sitk.LabelContour(reference_segmentation)

【SimpleITK】医疗影像分割结果评价指标计算

  1. sitk.SignedMaurerDistanceMap
reference_distance_map = sitk.Abs(sitk.SignedMaurerDistanceMap(reference_segmentation, squaredDistance=False, useImageSpacing=True))

【SimpleITK】医疗影像分割结果评价指标计算

  1. sitk.StatisticsImageFilter
statistics_image_filter = sitk.StatisticsImageFilter()
# Get the number of pixels in the reference surface by counting all pixels that are 1.
statistics_image_filter.Execute(reference_surface)
num_reference_surface_pixels = int(statistics_image_filter.GetSum()) 
print(num_reference_surface_pixels)

558

  1. sitk.LabelOverlapMeasuresImageFilter
overlap_measures_filter = sitk.LabelOverlapMeasuresImageFilter()
  1. sitk.HausdorffDistanceImageFilter
hausdorff_distance_filter = sitk.HausdorffDistanceImageFilter()