【SimpleITK】医疗影像分割结果评价指标计算_m0_37477175的博客-CSDN博客
chestnut-- 2020-01-16 17:20:48
文章目录
Overlap Measures
Overlap Measures
jaccard
dice
volume_similarity
false_negative
false_positive
其中:
Volume measures:
volumeSimilarityv1+v22∗(v1−v2)
建立枚举对象:
from enum import Enum# Use enumerations to represent the various evaluation measuresclass OverlapMeasures(Enum): jaccard, dice, volume_similarity, false_negative, false_positive = range(5) reference_segmentation = reference_segmentation_STAPLE12345
建立空的数组,为了后面保存结果:
# Empty numpy arrays to hold the results overlap_results = np.zeros((len(segmentations),len(OverlapMeasures.__members__.items()))) overlap_results123
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_results12345678910
结果:
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_STAPLE123
建立空的数组,为了后面保存结果:
surface_distance_results = np.zeros((len(segmentations),len(SurfaceDistanceMeasures.__members__.items()))) surface_distance_results12
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 = 1reference_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()) 12345678910
计算指标:
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)12345678910111213141516171819202122232425262728293031323334353637
结果:
[[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]]
为了理解上面的代码,我们需要了解豪斯多夫距离(hausdorff distance):
豪斯多夫距离
假设有两组集合A={a1,…,ap},B={b1,…,bq},A点集有p个点,而B点集有q个点。则这两个点集合之间的Hausdorff距离定义为:
H(A,B)=max(h(A,B),h(B,A))
其中,h(A,B)为,在A中的每一个点ai,到距离此点最近的B集合中的点的距离,然后对这些距离∣ai−bi∣进行排序。
A中有p个点,每个点都有B中距离最近的点。两两之间求距离,需要计算 p∗q次。
取所有距离中的最大值为h(A,B)的值。
some API
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)1
sitk.SignedMaurerDistanceMap
reference_distance_map = sitk.Abs(sitk.SignedMaurerDistanceMap(reference_segmentation, squaredDistance=False, useImageSpacing=True))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)12345
558
sitk.LabelOverlapMeasuresImageFilter
overlap_measures_filter = sitk.LabelOverlapMeasuresImageFilter()1
sitk.HausdorffDistanceImageFilter
hausdorff_distance_filter = sitk.HausdorffDistanceImageFilter()