Compute symmetric surface distances [Python]

Question

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))

Does this answer your question? [Evaluate neural network using predicted image and ground truth image](https://stackoverflow.com/questions/61575961/evaluate-neural-network-using-predicted-image-and-ground-truth-image) — Jingnan Jia, Sep 21 '22 at 21:13

score 4 · Accepted Answer · answered Nov 14 '17 at 20:25

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

do I have to recompute the Mauerer distances again, this time from the surface of the mask to the reference segmentation? then execute again the `label_intensity_statistics_filter` between contour of the reference and the Mauere distances of the mask? — RMS, Nov 15 '17 at 09:26
yes, you are asumming right. I am using that SimpleITK Notebook — RMS, Nov 15 '17 at 09:29
The need for symmetry was also raised by [this question](https://discourse.itk.org/t/root-mean-square-surface-distance-mauerer-distances-map-how-to-calculate-it/424) on the discourse forum. If this is not you (under a different name), the answer includes a link to updated code for symmetric surface distances. — zivy, Nov 16 '17 at 17:01

score 2 · Answer 2 · answered Nov 14 '17 at 13:28

Signed distance map is not symmetric. Hausdorff distance should be.

For mesh comparisons, I used metro in the past.

For Maurer, positive distances mean outside and negative distances mean inside. You should take absolute value if you want to calculate disagreement.

score 0 · Answer 3 · answered Sep 21 '22 at 14:32

seg-metrics is a simple package to compute different metrics for Medical image segmentation (images with suffix .mhd, .mha, .nii, .nii.gz or .nrrd), and write them to csv file.

Apaft from mean surface distance (MSD), there are some surface distance based metrics.

Hausdorff distance
Hausdorff distance 95% percentile
Mean (Average) surface distance
Median surface distance
Standard deviation of surface distance

Note: These metrics are symmetric, which means the distance from A to B is the same as the distance from B to A.

https://pypi.org/project/seg-metrics/

An example:

Install this package at first. pip install seg-metrics

import seg_metrics.seg_metrics as sg

labels = [0, 4, 5 ,6 ,7 , 8]
gdth_file = 'data/gdth.mhd'  # ground truth image full path
pred_file = 'data/pred.mhd'  # prediction image full path
csv_file = 'metrics.csv'

metrics = sg.write_metrics(labels=labels[1:],  # exclude background
                  gdth_path=gdth_file,
                  pred_path=pred_file,
                  csv_file=csv_file, 
                  metrics=['hd', 'msd'])

Compute symmetric surface distances [Python]

3 Answers3