Computing 95% Hausdorff Distance

Beginner Questions
,
1

Hello,

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: 34_Segmentation_Evaluation . 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!

2

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
3

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.

cheers

import numpy as np
import SimpleITK as sitk

NUM_LABELS = 18

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
4

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

Hi @zivy

It seems the sample code you’ve provided calculates the boundary distances, is there an equivalent way to calculate the point-based differences so that we can calculate the 95% HD based on these?

Also if you could clear up my understanding of the following that would be appreciated. From my empirical tests, the sitk.DirectedHausdorffDistanceImageFilter calculates the HD using the point-wise distances, not the boundary distances. However, in the ITK docs, it’s stated that this filter makes use of the sitk.SignedMaurerDistanceMap (ref: https://itk.org/Doxygen/html/classitk_1_1DirectedHausdorffDistanceImageFilter.html). As I understand it, the SignedMaurerDistanceMap calculates distances to the boundary, so how can this be used for the point-wise calc?

Thanks!

6

If you have a look at the code, you can see that the filter takes the max( signed_distance, 0), i.e. distances inside the object region are zero.

See line 174:

// do the work
  while (!it1.IsAtEnd())
  {
    if (Math::NotExactlyEquals(it1.Get(), NumericTraits<InputImage1PixelType>::ZeroValue()))
    {
      // The signed distance map is calculated, but we want the calculation based on the
      // unsigned distance map.  Therefore, we set all distance map values less than 0 to 0.
      const RealType val2 = std::max(static_cast<RealType>(it2.Get()), NumericTraits<RealType>::ZeroValue());
      maxDistance = std::max(maxDistance, val2);
      sum += val2;
      ++pixelCount;
    }
7

I guess following code should work to compute the symmetric point-wise hausdorff distance (using itk pythonic api):

def hausdorff_pointwise_distance(y_pred: itkImage, y_ref: itkImage) -> Dict[str, float]:
    """Compute symmetric point-wise distances between two binary masks

    Args:
        y_pred (itkImage): predicted segmentation
        y_ref (itkImage): reference segmentation

    Returns:
        Dict[str, float]: keys are 'mean', 'median', 'std', 'max'
    """
    seg_distance = itk.signed_maurer_distance_map_image_filter(
        y_pred, use_image_spacing=True, squared_distance=False, inside_is_positive=False
    )
    ref_distance = itk.signed_maurer_distance_map_image_filter(
        y_ref, use_image_spacing=True, squared_distance=False, inside_is_positive=False
    )

    # get distance inside foreground label
    ref2seg_distance = seg_distance[np.nonzero(y_ref)]
    seg2ref_distance = ref_distance[np.nonzero(y_pred)]

    # compute statistics on symmetric distances (both directions)
    all_surface_distances = np.concatenate(
        (ref2seg_distance, seg2ref_distance), axis=None
    )
    all_surface_distances[all_surface_distances <= 0.0] = 0.0

    mean_surface_distance = np.mean(all_surface_distances)
    median_surface_distance = np.median(all_surface_distances)
    std_surface_distance = np.std(all_surface_distances)
    max_surface_distance = np.max(all_surface_distances)
    return {
        "mean": mean_surface_distance,
        "median": median_surface_distance,
        "std": std_surface_distance,
        "max": max_surface_distance,
    }