Compute symmetric surface distances [Python]

I want to compute surface distance metrics between 2 binary objects, aka segmentations of liver tumors. I am looking to compute:

Average symmetric surface distance
Root mean square symmetric distance
Hausdorff distance (also known as maximum symmetric distance)

I found two libraries that could help me compute those metrics, but I am getting conflicting results, so i am confused how they work.

This is my code for Simple ITK and MedPy.

from medpy import metric
import pandas as pd
import SimpleITK as sitk
import numpy as np
reference_segmentation = sitk.ReadImage('tumorSegm', sitk.sitkUInt8)
segmentation = sitk.ReadImage('tumorSegm2',sitk.sitkUInt8)
class SurfaceDistanceMeasuresITK(Enum):
    hausdorff_distance, max_surface_distance, avg_surface_distance, median_surface_distance, std_surface_distance = range(5)

class MedpyMetricDists(Enum):
    hausdorff_distance, avg_surface_distance, avg_symmetric_surface_distance = range(3)


  surface_distance_results = np.zeros((1,len(SurfaceDistanceMeasuresITK.__members__.items())))
surface_dists_Medpy = np.zeros((1,len(MedpyMetricDists.__members__.items())))
segmented_surface = sitk.LabelContour(segmentation)

# init signed mauerer distance as reference metrics
reference_distance_map = sitk.Abs(sitk.SignedMaurerDistanceMap(reference_segmentation, squaredDistance=False, useImageSpacing=True))

label_intensity_statistics_filter = sitk.LabelIntensityStatisticsImageFilter()
label_intensity_statistics_filter.Execute(segmented_surface, reference_distance_map)

hausdorff_distance_filter = sitk.HausdorffDistanceImageFilter()
hausdorff_distance_filter.Execute(reference_segmentation, segmentation)

surface_distance_results[0,SurfaceDistanceMeasuresITK.hausdorff_distance.value] = hausdorff_distance_filter.GetHausdorffDistance()
surface_distance_results[0,SurfaceDistanceMeasuresITK.max_surface_distance.value] = label_intensity_statistics_filter.GetMaximum(label)
surface_distance_results[0,SurfaceDistanceMeasuresITK.avg_surface_distance.value] = label_intensity_statistics_filter.GetMean(label)
surface_distance_results[0,SurfaceDistanceMeasuresITK.median_surface_distance.value] = label_intensity_statistics_filter.GetMedian(label)
surface_distance_results[0,SurfaceDistanceMeasuresITK.std_surface_distance.value] = label_intensity_statistics_filter.GetStandardDeviation(label)

surface_distance_results_df = pd.DataFrame(data=surface_distance_results, index = list(range(1)),
                              columns=[name for name, _ in SurfaceDistanceMeasuresITK.__members__.items()])

img_array = sitk.GetArrayFromImage(reference_segmentation)
seg_array = sitk.GetArrayFromImage(segmentation)
# reverse array in the order x, y, z
img_array_rev = np.flip(img_array,2)
seg_array_rev = np.flip(seg_array,2)
vxlspacing = segmentation.GetSpacing()

surface_dists_Medpy[0,MedpyMetricDists.hausdorff_distance.value] = metric.binary.hd(seg_array_rev,img_array_rev, voxelspacing=vxlspacing)
surface_dists_Medpy[0,MedpyMetricDists.avg_surface_distance.value] = metric.binary.asd(seg_array_rev,img_array_rev, voxelspacing=vxlspacing)
surface_dists_Medpy[0,MedpyMetricDists.avg_symmetric_surface_distance.value] = metric.binary.assd(seg_array_rev,img_array_rev, voxelspacing=vxlspacing)

surface_dists_Medpy_df = pd.DataFrame(data=surface_dists_Medpy, index = list(range(1)),
                              columns=[name for name, _ in MedpyMetricDists.__members__.items()])

At first sight I don't think that SimpleITK computes the symmetric distances. Is there any implementation for those in that library? How can I obtain them?
Is MedPy a reliable library? Can I calculate the symmetric root mean square with it?
Other recommendations of libraries for computing surface distance metrics?
Should I compute the Absolute Values for the Mauerer Distances Map? I am not sure how it affects the results. reference_distance_map = sitk.Abs(sitk.SignedMaurerDistanceMap(reference_segmentation, squaredDistance=False, useImageSpacing=True))

Solution

@Roxanne

I am assuming here that you are being confused by the surface distance measures computed in this SimpleITK notebook?

The rest of the answer refers to that code.

The mean/std/median/maximum are not symmetric (Hausdorff is).

Using SimpleITK you can compute the symmetric mean and standard deviation by computing the mean and standard deviation for the segmentation and then for the reference (the code does it for the segmentation, so just switch roles and you get it for the reference).

Now you have the mean and standard deviations from two samples. To get the size of a sample just call:

label_intensity_statistics_filter.GetNumberOfPixels(label)

To compute the symmetric mean and standard deviation from knowlege of n1,m1,s1, n2,m2,s2:

m = (n1*m1 + n2*m2)/(n1+n2)
s = np.sqrt((n1*(s1**2+(m1-m)**2) + n2*(s2**2+(m2-m)**2))/(n1+n2))

Note that the sample estimate for the standard deviation is the biased version (similar to the default behavior of numpy.std).

If you have additional questions please post to the ITK discourse forum.