Calculating mutual information of two images

Beginner Questions
,
1

I have corresponding MRI and CT image volumes. I have registered them using the mattes mutual information metric. The results seem fine but I would like to be able to calculate the metric of two 3D images that are already registered. As of now the only way to calculate the metric I know of is to call the MetricEvaluate method of the ImageRegistrationMethod class. I have written a function to calculate the mutual information given the two image volumes (before resampling) and the transformation found during registration:

def mutual_info(transform, mr, ct):
    registration_method = sitk.ImageRegistrationMethod()
    registration_method.SetMetricAsMattesMutualInformation(numberOfHistogramBins=50)
    registration_method.SetMetricSamplingStrategy(registration_method.RANDOM)
    registration_method.SetMetricSamplingPercentage(0.01)
    registration_method.SetInterpolator(sitk.sitkLinear)
    registration_method.SetInitialTransform(transform)
    return registration_method.MetricEvaluate(mr, ct)

But there are two problems:

  1. The resulting value is different from the one achieved at the end of the registration process. Why is this the case?
  2. I would like to calculate the mutual information of two images that were already resampled according to the transformation found during registration. And I would also like to calculate the mutual information of only the upper part of the images. Is there a way to simply calculate the mutual information of two images without a transformation?
2

Hello @Flo,

  1. The value is different because of the use of random sampling. The way the function is written, you will obtain a different value every time it is run (additional details about registration variability on Read-The-Docs).
  2. When you say “no transformation” you are actually saying “the identity transformation”, so we always have a transformation, explicit or implicit.

Recommendations:
a. Sampling strategy of NONE.
b. Global number of threads - one.
c. Metric value in a sub image, first 20 slices of CT, assuming ct_image is the fixed image, mutual_info(transform, ct_image[:,:,0:20], mr_image) should work.

2 Likes
3

Thanks a lot for your answer!
I tried to implement your suggestions. In the beginning of my jupiter notebook I run this line:

sitk.ProcessObject.SetGlobalDefaultNumberOfThreads(1)

My mutual_info function now looks like this:

def mutual_info(transform, mr, ct):
    registration_method = sitk.ImageRegistrationMethod()
    registration_method.SetMetricAsMattesMutualInformation(numberOfHistogramBins=50)
    registration_method.SetMetricSamplingStrategy(registration_method.NONE)
    registration_method.SetInterpolator(sitk.sitkLinear)
    registration_method.SetInitialTransform(transform)
    return registration_method.MetricEvaluate(mr, ct)

But this still results in different values than than using registration_method.GetMetricValue() directly after registration. Here Is the code I use to register the images:

def register(fixed, moving, multi_res=False, initial_transform=None):

    registration_method = sitk.ImageRegistrationMethod()
    registration_method.SetMetricAsMattesMutualInformation(numberOfHistogramBins=50)
    registration_method.SetMetricSamplingStrategy(registration_method.RANDOM)
    registration_method.SetMetricSamplingPercentage(0.01)
    registration_method.SetInterpolator(sitk.sitkLinear)

    # Optimizer settings.
    registration_method.SetOptimizerAsGradientDescent(learningRate=1.0, numberOfIterations=1000, convergenceMinimumValue=1e-6, convergenceWindowSize=20) #estimateLearningRate=sitk.ImageRegistrationMethod.Never)
    registration_method.SetOptimizerScalesFromPhysicalShift()

    # Setup for the multi-resolution framework.
    if multi_res:
        registration_method.SetShrinkFactorsPerLevel(shrinkFactors=[4, 2, 1])
        registration_method.SetSmoothingSigmasPerLevel(smoothingSigmas=[2, 1, 0])
        registration_method.SmoothingSigmasAreSpecifiedInPhysicalUnitsOn()

    if initial_transform is None:
        initial_transform = sitk.CenteredTransformInitializer(fixed,
                                                          moving,
                                                          sitk.Euler3DTransform(),
                                                          sitk.CenteredTransformInitializerFilter.GEOMETRY)

    registration_method.SetInitialTransform(initial_transform)
    final_transform = registration_method.Execute(fixed, moving)

    print('Final metric value: {0}'.format(registration_method.GetMetricValue()))
    print('Optimizer\'s stopping condition, {0}'.format(registration_method.GetOptimizerStopConditionDescription()))

    registration_method.SetMetricSamplingStrategy(registration_method.NONE)
    metric_1 = registration_method.GetMetricValue()
    metric_2 = mutual_info(final_transform, fixed, moving)

    moving_resampled = sitk.Resample(moving, fixed, final_transform, sitk.sitkLinear, 0.0, moving.GetPixelID())
    identity_transform = sitk.CenteredTransformInitializer(fixed,
                                                           moving_resampled,
                                                           sitk.Euler3DTransform(),
                                                           sitk.CenteredTransformInitializerFilter.GEOMETRY)

    metric_3 = mutual_info(identity_transform, fixed, moving_resampled)

    return final_transform, metric_1, metric_2, metric_3

After Running this code I get:

metric_1 = -0.5781037066906896
metric_2 = -0.5485657392397647
metric_3 = -0.5086132295711404

If I run the code repeatedly metric_1 is always somewhere between -0.57 and -0.60 while metric_2 is always close to -0.55

Why are they still different?

I added metric_3 to show the case where I use my ‘mutual_info’ function to calculate the metric after the CT was resampled by creating an identity transformation. If everything was correct then I would expect to find metric_1 == metric_2 == metric_3. But I always get metric_1 > metric_2 > metric_3.

Do you have any further suggestions?

4

Hello @Flo,

Each of these statements is estimating the metric with a different point set:

metric_1: The last metric value computed during registration using the last sampling. Simply setting the sampling strategy after the registration was completed to None does not cause it to compute the metric again. The GetMetricValue just returns the value, no computation is performed (same value as printed in the “Final metric value:” print statment).

metric_2: The metric is computed using only fixed image points that map inside the moving image, those that fall outside are ignored.

metric 3: Because the moving image was resampled, points that were originally mapped outside the moving image are now inside the moving_resampled image and have a value of 0.0.

One last thing, to create the identity transformation, just do identity_transform = sitk.Transform().

2 Likes
5

Thank you so much @zivy

I understand now.