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?