title: Compute symmetric surface distances [Python]
body: <p>I want to compute surface distance metrics between 2 binary objects, aka segmentations of liver tumors. I am looking to compute:</p>
  
  <ul>
  <li>Average symmetric surface distance</li>
  <li>Root mean square symmetric distance</li>
  <li>Hausdorff distance (also known as maximum symmetric distance)</li>
  </ul>
  
  <p>I found two libraries that could help me compute those metrics, but I am getting conflicting results, so i am confused how they work.</p>
  
  <ul>
  <li>SimpleITK:<a href="https://itk.org/SimpleITKDoxygen/html/classitk_1_1simple_1_1SignedMaurerDistanceMapImageFilter.html" rel="nofollow noreferrer">https://itk.org/SimpleITKDoxygen/html/classitk_1_1simple_1_1SignedMaurerDistanceMapImageFilter.html</a>
  <a href="https://itk.org/SimpleITKDoxygen/html/classitk_1_1simple_1_1HausdorffDistanceImageFilter.html" rel="nofollow noreferrer">https://itk.org/SimpleITKDoxygen/html/classitk_1_1simple_1_1HausdorffDistanceImageFilter.html</a></li>
  <li>MedPy:
  <a href="https://github.com/loli/medpy/" rel="nofollow noreferrer">https://github.com/loli/medpy/</a>. <a href="http://loli.github.io/medpy/metric.html" rel="nofollow noreferrer">DOCS</a></li>
  </ul>
  
  <p>This is my code for Simple ITK and MedPy.</p>
  
  <pre><code>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()])
  </code></pre>
  
  <ol>
  <li><p>At first sight I don't think that SimpleITK computes the
  <strong>symmetric</strong> distances. Is there any implementation for those in that library? How can I obtain them?</p></li>
  <li><p>Is MedPy a reliable library? Can I calculate the symmetric root mean
  square with it?</p></li>
  <li>Other recommendations of libraries for computing surface distance
  metrics?</li>
  <li>Should I compute the Absolute Values for the Mauerer Distances Map? I am not sure how it affects the results.   <code>reference_distance_map = sitk.Abs(sitk.SignedMaurerDistanceMap(reference_segmentation, squaredDistance=False, useImageSpacing=True))
  </code></li>
  </ol>