Dice score and Hausdorff distance for inter-modal image registration

I have 3D CT and MR images of the brain. I want to register these images and find out Dice score, Hausdorff distance, Jaccard index as measure of registration accuracy.
I have segmented the whole brain (both CT and MR images). I have also prepared binary masks of the whole brain (both CT and MR). Samples of segmentations and masks are shown below.

image image

image image

I studied a similar registration code given in SimpleITK notebooks/65_Registration_FFD.ipynb. This code is for intra-modal lung registration. It uses a utility called lung_label, to extract the lung from the whole mask.

How can I prepare a similar label for my case?
Is there any example code for inter-modal registration, where Dice score and Hausdorff distance are computed?

I am using my own metric for registration. So I cannot use an existing registration algorithm for the job.

Hello @debapriya,

Not sure what is the problem you are encountering. The computation of the segmentation measures is independent of the inter/intra modal aspect of the registration or the type of transformation estimated via registration. All that is required is the estimated transformation.

The pseudo-code for what you want to do is:

fixed_image = 
fixed_image_segmentation = 
moving_image = 
moving_image_segmentation = 
# Whatever the segmentation label of interest is. There can be multiple labels/structures in the segmentation image so identify which one we are interested in. 
# Value is likely 1 but can be any arbitrary value e.g. 42.
segmentation_label_of_interest =  

# register the images and resample the moving_image_segmentation
tx = register(fixed_image, moving_image)
moving_image_segmentation_pre_registration = sitk.Resample(moving_image_segmentation,
                                                    fixed_image_segmentation,
                                                    sitk.Transform(),
                                                    sitk.sitkNearestNeighbor)

moving_image_segmentation_post_registration = sitk.Resample(moving_image_segmentation,
                                                    fixed_image_segmentation,
                                                    tx,
                                                    sitk.sitkNearestNeighbor)

label_overlap_measures_filter = sitk.LabelOverlapMeasuresImageFilter()
label_overlap_measures_filter.Execute(fixed_image_segmentation, moving_image_segmentation_pre_registration)
print(
f"Dice coefficient before registration: {label_overlap_measures_filter.GetDiceCoefficient(segmentation_label_of_interest):.2f}"
)
label_overlap_measures_filter.Execute(fixed_image_segmentation, moving_image_segmentation_post_registration)
print(
f"Dice coefficient after registration: {label_overlap_measures_filter.GetDiceCoefficient(segmentation_label_of_interest):.2f}"
)
...

2 Likes

Thank you @zivy for this lucid explanation. It is very helpful.

Following your instructions, I set segmentation_label_of_interest = 100 and executed the code. I am running into the following error.

My fixed_image size is 512, 512, 26 and moving_image size is 276, 384, 21. fixed_mask and moving_mask sizes are same as fixed_image and moving_image sizes respectively.

Hello @debapriya,

There was an assumption in the original pseudo-code that the fixed and moving segmentations had the same metadata (origin, spacing…). This was incorrect.

We need to resample the original moving image segmentation onto the fixed image grid using the identity transform so that we can compare the fixed and moving segmentations prior to registration. The pseudo-code was updated accordingly. Please see the updated code above.

Thank you @zivy . I updated the code. Now I am getting a different error about pixel type not being supported.

Code for loading data and masks is given below

#Load image

fixed_image = sitk.ReadImage("Data/M_CT.mha", sitk.sitkFloat64)
moving_image = sitk.ReadImage("Data/M_MR.mha", sitk.sitkFloat64)

# Load mask

fixed_mask = sitk.ReadImage("Data/M_CT_Mask.mha", sitk.sitkFloat64)
moving_mask = sitk.ReadImage("Data/M_MR_Mask.mha", sitk.sitkFloat64)

I tried with sitkFloat32 as well. I am getting the same error.

Hello @debapriya,

A label image is an image with integer types (sitkUInt8, sitkUInt16…). The way the masks are read in the provided code snippet is incorrect. The easiest is to not specify the pixel type and if it is indeed a label image it will be one of the integer types:

fixed_mask = sitk.ReadImage("Data/M_CT_Mask.mha")
moving_mask = sitk.ReadImage("Data/M_MR_Mask.mha")
1 Like

How can I make sure that my mask is a label image?

I am using Slicer 5.4.0 to create the mask. I am selecting 0 as Outside fill value and 100 as Inside fill value. I am choosing Create New LabelMapVolume as... from the dropdown.

I am saving the data as a .mha volume

image

I am still getting the following error

Hello @debapriya,

I suspect you are saving it as a labelmap which is different from a label image. Pixel type of sitkLabelUInt8 vs. sitkUInt8 (see all pixel types here). Try casting to sitkUInt8 when reading (e.g. fixed_mask = sitk.ReadImage("Data/M_CT_Mask.mha", sitk.sitkUInt8)).

If this doesn’t address the issue, please share a sample dataset so that we can help you resolve the issue.

Casting to sitkUInt8 solved my purpose. Thank you @zivy.