Overlap Measures

• jaccard
• dice
• volume_similarity
• false_negative
• false_positive

其中：
Volume measures:
$$volumeSimilarity\frac{2\ast (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]]

为了理解上面的代码，我们需要了解豪斯多夫距离（hausdorff distance）：

豪斯多夫距离

参考维基百科, 博客。

假设有两组集合 A = { a 1 , … , a p } A=\{a_1,…,a_p\} A={a1​,…,ap​}， B = { b 1 , … , b q } B=\{b_1,…,b_q\} B={b1​,…,bq​}，A点集有p个点，而B点集有q个点。则这两个点集合之间的Hausdorff距离定义为:
$$H(A,B)=max(h(A,B),h(B,A))$$
其中，$$h(A,B)$$为，在A中的每一个点$a_i$，到距离此点最近的B集合中的点的距离，然后对这些距离$|a_i-b_i|$进行排序。
A中有p个点，每个点都有B中距离最近的点。两两之间求距离，需要计算 $p\ast q$次。
取所有距离中的最大值为$h(A,B)$的值。

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)

2. sitk.SignedMaurerDistanceMap

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

3. 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

4. sitk.LabelOverlapMeasuresImageFilter

overlap_measures_filter = sitk.LabelOverlapMeasuresImageFilter()

5. sitk.HausdorffDistanceImageFilter

hausdorff_distance_filter = sitk.HausdorffDistanceImageFilter()