# Calculating mutual information of two images

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?

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

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

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

Thank you so much @zivy

I understand now.