Computing 95% Hausdorff Distance


I am using Python and SimpleITK, and I would like to access the raw distance values which are used in the calculation of the Hausdorff Distance in order to calculate the 95% Hausdorff Distance. When I use sitk.HausdorffDistanceImageFilter I am able to retrieve the maximum distance or the average distance. After looking around at the itk documentation, what I need might be located in itk’s DirectedHausdorffDistanceImageFilter. However, I don’t know how to access that output from within Python using SimpleITK.

I found this documentation here: . However, it shows how to get the raw distance values for a surface-based Hausdorff distance, while the sitk filter uses point based segmentation results. It mentions how these two methods of calculation are not equivalent:

“the Hausdorff distance for the contour/surface representation and the discrete point set representing the segmented object differ, and that there is no correlation between the two.”

After looking more into the documentation, it seems like I could write something in Python using sitk.SignedMaurerDistanceMap, and looping through the voxels to find the maximum, but it seems kind of inefficient if something optimized in the c++ libraries already exists.

Is there a way to access the “raw” distance values so that I can compute the 95% Hausdorff Distance (or any other percentile, for that matter) and leverage the speed of the itk libraries and use a discrete point set segmented object?

Thank you in advance for your help!

Hello @emckenzi123,

There is no way to access the raw distance values from the HausdorffDistanceImageFilter. Also, the Hausdorff distance definition is a single number (max of the two non-symmetric distances, unclear what you mean with respect to percentile in that context).

If you want to compute various statistics on the distances between the two point sets that is readily done via a distance map. Add the following as a last code cell to the notebook you refered to. The code computes the Hausdorff distance and you have access to the discrete set of distances used for the computation:

seg_distance_map = sitk.Abs(sitk.SignedMaurerDistanceMap(seg, squaredDistance=False, useImageSpacing=True))
reference_segmentation1_distance_map = sitk.Abs(sitk.SignedMaurerDistanceMap(reference_segmentation1, squaredDistance=False, useImageSpacing=True))
reference_segmentation2_distance_map = sitk.Abs(sitk.SignedMaurerDistanceMap(reference_segmentation2, squaredDistance=False, useImageSpacing=True))

dist_1_seg = sitk.GetArrayViewFromImage(seg_distance_map)[sitk.GetArrayViewFromImage(reference_segmentation1)==1]
dist_seg_1 = sitk.GetArrayViewFromImage(reference_segmentation1_distance_map)[sitk.GetArrayViewFromImage(seg)==1]
print(f'Hausdorff distance: {np.concatenate([dist_1_seg, dist_seg_1]).max()}')

dist_2_seg = sitk.GetArrayViewFromImage(seg_distance_map)[sitk.GetArrayViewFromImage(reference_segmentation2)==1]
dist_seg_2 = sitk.GetArrayViewFromImage(reference_segmentation2_distance_map)[sitk.GetArrayViewFromImage(seg)==1]
np.concatenate([dist_2_seg, dist_seg_2]).max()
print(f'Hausdorff distance: {np.concatenate([dist_2_seg, dist_seg_2]).max()}')
1 Like

Thank you zivy for your help. I was trying to use sitk to calculate the 95% Hausdorff Distance as defined by the Plastimatch software (see Plastimatch Hausdorff Documentation ). I was not able to recreate their values using sitk, however when digging deeper into their source code, I discovered that their default distance algorithm was a Signed Danielsson Distance Map. There were a few other implementation details that I also discovered accounted for my discrepancies. Anyways, for future people trying to recreate the plastimatch 95% Hausdorff Distance (boundary) in sitk, I’ve included my code below. I ended up using Maurer Distance for speed, but if you want an exact match, I left the code for using the Danielsson distance commented.


import numpy as np
import SimpleITK as sitk


in_img = sitk.ReadImage('seg_contour.nii.gz')
in_img = sitk.GetArrayFromImage(in_img)
ref_img = sitk.ReadImage('ref_contour.nii.gz')
ref_img = sitk.GetArrayFromImage(ref_img)

for label in range(1,NUM_LABELS):
    ### Get binary mask of the current label
    in1 = np.zeros(in_img.shape)
    ref1 = np.zeros(ref_img.shape)
    in1[in_img==label] = 1.0
    ref1[ref_img==label] = 1.0
    ref_contour_img = sitk.GetImageFromArray(ref1)
    test_contour_img = sitk.GetImageFromArray(in1)
    reference_segmentation = sitk.Cast(ref_contour_img, sitk.sitkUInt32)
    seg = sitk.Cast(test_contour_img, sitk.sitkUInt32)
    reference_surface = sitk.LabelContour(reference_segmentation, False)
    seg_surface = sitk.LabelContour(seg, False)

    ### Get distance map for contours (the distance map computes the minimum distances) 

    # seg_distance_map = sitk.Abs(sitk.SignedDanielssonDistanceMap(seg_surface, insideIsPositive=False, squaredDistance=False, useImageSpacing=True))
    seg_distance_map = sitk.Abs(sitk.SignedMaurerDistanceMap(seg_surface, squaredDistance=False, useImageSpacing=True))
    # reference_segmentation_distance_map = sitk.Abs(sitk.SignedDanielssonDistanceMap(reference_surface, insideIsPositive=False, squaredDistance=False, useImageSpacing=True))
    reference_segmentation_distance_map = sitk.Abs(sitk.SignedMaurerDistanceMap(reference_segmentation, squaredDistance=False, useImageSpacing=True))

    ### Find the distances to surface points of the contour.  Calculate in both directions         
    dist_seg = sitk.GetArrayViewFromImage(seg_distance_map)[sitk.GetArrayViewFromImage(reference_surface)==1]
    dist_ref = sitk.GetArrayViewFromImage(reference_segmentation_distance_map)[sitk.GetArrayViewFromImage(seg_surface)==1]

    ### Find the 95% Distance for each direction and average        
    print((np.percentile(dist_ref, 95) + np.percentile(dist_seg, 95)) / 2.0)
1 Like

Nice work!

I might be mistaken but you are not using the image spacing information.
Please review the changes I made:

from functools import partial

import numpy as np
import SimpleITK as sitk
from SimpleITK import GetArrayViewFromImage as ArrayView

prediction = sitk.ReadImage("/path/to/prediction")
gold = sitk.ReadImage("/path/to/gold")

distance_map = partial(sitk.SignedMaurerDistanceMap, squaredDistance=False, useImageSpacing=True)

num_labels = 18
for label in range(1, num_labels + 1):
    gold_surface = sitk.LabelContour(gold == label, False)
    prediction_surface = sitk.LabelContour(prediction == label, False)

    ### Get distance map for contours (the distance map computes the minimum distances)
    prediction_distance_map = sitk.Abs(distance_map(prediction_surface))
    gold_distance_map = sitk.Abs(distance_map(gold_surface))

    ### Find the distances to surface points of the contour.  Calculate in both directions
    gold_to_prediction = ArrayView(prediction_distance_map)[ArrayView(gold_surface) == 1]
    prediction_to_gold = ArrayView(gold_distance_map)[ArrayView(prediction_surface) == 1]

    ### Find the 95% Distance for each direction and average
    print((np.percentile(prediction_to_gold, 95) + np.percentile(gold_to_prediction, 95)) / 2.0)