Inconsistent Mattes Mutual Information metric Value.

Algorithms
,
1

I am trying to register MRI- CT volumes from " Report of the AAPM Radiation Therapy Committee Task Group No. 132" (dataset). The dataset has CT, MR and PET volumes of geometric and anatomical phantoms. The offsets along all axis between the volumes are precisely known. My registration algorithm is converging and results are very close to ground truth but needs some more tuning(Optimizer Regular step gradient descent with Mattes Mutual Information as metric). On close examination I observed that the metric values with ground truth transformation is less than results with my algorithm. It should not be the case. Following example will clarify my problem.

Registration : Phantom Dataset 2 CT volume with Phantom Dataset 1 MRI T2 volume (Rigid registration with translation only) Link

CT Volume size: 512 512 46
CT Volume spacing: 0.7031487 0.7031487 3
CT Volume Direction: 1 0 0 0
CT Volume Origin: -180.23 -180.23 -66.5
MR T2 Volume size: 512 512 46
MR T2 Volume spacing: 0.7031487 0.7031487 3
MR T2 Volume Direction: 1 0 0 0
MR T2 Volume Origin: -180.23 -180.23 -66.5

Ground truth: Offsets in CT, left = 1.0 cm, to anterior = 0.5 cm, to inferior = 1.5 cm
Ground truth Metric Value: 0.8699 (bin size 100)
My registration results:
Center: [-0.575507, -0.575507, 1]
Translation: [-10.2376, 5.28658, -15.0103]
Metric Value: 0.8819 (bin size 100)

Following is how I am estimating the Metric value

Euler3DTransform regTransform = new Euler3DTransform();
regTransform.SetParameters(new VectorDouble() { 0, 0, 0, -10.0, 5.0, -15.0 });
regTransform.SetCenter(new VectorDouble() { -0.575507, -0.575507, 1 });
uint numberOfBins = 100

ImageRegistrationMethod R = new ImageRegistrationMethod();
R.SetMetricAsMattesMutualInformation(numberOfBins);
R.SetMetricSamplingStrategy(ImageRegistrationMethod.MetricSamplingStrategyType.NONE);
R.SetInterpolator(InterpolatorEnum.sitkLinear);
R.SetInitialTransform(regTransform);
double metricValue = R.MetricEvaluate(fixedImage, movingImage);

Even with other metric i.e. Joint Histogram Mutual Information I am facing same issue.
6.2789 (Ground truth) vs 6.3659 (My registration)

Am I missing something here? Could this be the reason registration is not converging towards ground truth?

2

@zivy Could you please look at this query.
Thanks in advance.

3

Hello @gaurav_patel,

What you are seeing is common in optimization of non-convex functions. Registration is just a specific use case of this general class of problems. The optimization converges to a local extreme point of the optimized function, in ITK/SimpleITK a local minimum. This is why registration success is so dependent on initialization.

In your particular case, the results are very close to the known ground truth in parameter space. This is for a single registration run. If you run the registration multiple times, even when you start with the same parameters you will converge to slightly different results due to the stochastic nature of the process (read more here).

There is no way to ensure that you have ended up at the desired minimum point in the parameter space. Once you reach the local minimum you can explore the region in parameter space around it, similar to what you did. Set a grid in the parameter space and evaluate the optimized function for those parameter values (see Exhaustive optimizer in this notebook).

Finally, the errors in parameter space do not mean much. They are useful for technical analysis and look great in publications, but should be analyzed cautiously. The relevant quantity is Target Registration Error, distance between corresponding points in the two coordinate systems after one of them was mapped to the other coordinate system. The reason the errors in parameter space don’t mean much is that the same error in parameter space may have significantly different effects in physical space. The effect of a rotation error of 1 degree depends on the location of the center of rotation with respect to the target region. If the center of rotation is close to the target the effect of the error is smaller. In the extreme, if it coincides with the target, the rotational errors have no effect.

Finally, evaluating results is context dependent. The same results may be sufficiently accurate for task A, but not for task B.

See this book chapter for a more general discussion of registration evaluation.

2 Likes
4

Thank you for detailed response. I will look into all the resources you have shared.

There is no way to ensure that you have ended up at the desired minimum point in the parameter space. Once you reach the local minimum you can explore the region in parameter space around it, similar to what you did. Set a grid in the parameter space and evaluate the optimized function for those parameter values (see Exhaustive optimizer in this notebook).

I will try this out.

I have one specific issue that I am trying to figure out. I have mentioned in original query but I am reporting here again.

I am trying to register MRI- CT volumes from " Report of the AAPM Radiation Therapy Committee Task Group No. 132" (dataset ).

Registration : Phantom Dataset 2 CT volume with Phantom Dataset 1 MRI T2 volume (Rigid registration with translation only) Link

Specific problem I want your opinion on is following
Ground truth Translation [-10.0, 5.0, -15.0]
Ground truth Metric Value: 0.8699 (bin size 100)

For some other parameters
Translation: [-10.2376, 5.28658, -15.0103]
Metric Value: 0.8819 (bin size 100)

Translation: [-10.23, 5.0, -15.0]
Metric Value: 0.8736 (bin size 100)

Translation:[-10.0, 5.28, -15.0]
Metric Value: 0.8753 (bin size 100)

Translation:[-10.2, 5.2, -15.0]
Metric Value: 0.8772 (bin size 100)

How can metric value for some transformation parameters be more than ground truth transformation parameters? I have mentioned in my original query how I am utilizing simpleitk ImageRegistrationMethod object to calculate Mutual information metric.

5

I would recommend to create 2D surface plots of the similarity metric value around the known optimum transformation parameters, perturbin 2 transformation parameters each time (e.g., translation L, P; translation L, S; rotation LR, AP; rotation LR, IS). Typical issues that make the optimizer end up not finding the correct transformation:

  • very noisy surface with lots of local optima
  • peak has a large flat top
  • instead of a peak you have a long ridge

You need to try different metrics, metric parameters, sampling percentage, maybe image proprocessing, etc. - until you see metric surface plots that have only one peak, at the optimal transformation parameters.

If you already have one peak but the optimizer does not find it then you need to choose a different optimizer or optimization parameters or initial guess.

Note that most modern registration frameworks, such as Elastix or ANTs come with default presets that should work well for most cases. So, if you don’t have special, particularly challenging images then they should work without the need for tuning any parameters.

6

I had a look at the data sets and they are just some anatomical phantoms. The phantom contains a few major anatomical structures but they are made up of homogeneous materials. This kind of phantom is perfect for testing software workflows, visualization, application of transforms, etc.

image
image1471Ă—1128 144 KB

However, these images do not resemble real anatomical images in term of image content, details, and quality. You cannot use these data sets for evaluating automatic intensity-based image registration algorithms, because these algorithms rely on details that are present in real patient images but completely missing in these phantoms. Registration algorithms that work on real patient images may not work well on these phantom images; and registration algorithms that are developed to work on the phantom images may fail on real patient images.

Anyway, I gave it a try and run Elastix in 3D Slicer using default registration settings and I got visually perfect alignment:

This was the registration result:

    1  -0.000358915  -6.19205e-05  2.47632
    0.000358908  1  -0.000113049  -5.50905
    6.19611e-05  0.000113027  1  -12.1
    0  0  0  1

The computed transformation error is less than 0.5mm. Considering that the voxel sizes are 0.9x0.9x3.0 for the CT and 1.8x1.8x3.0 for the MRI, the difference of the registration compared to the ground truth was 1/4 of a voxel, which means that the registration was perfect.

I’ve even tried the default deformable (rigid+bspline) preset in Elastix and that worked, too. Deformable registration is expected to perform poorly in homogeneous regions, like those that you have in a phantom, but it still was OK. The displacement field only had ripples outside the phantom region (mostly where the “ImSimQA” annotation was burnt into the slices).

image
image1226Ă—1135 275 KB

1 Like
7

Thank you for such detailed response.

I would recommend to create 2D surface plots of the similarity metric value around the known optimum transformation parameters, perturbin 2 transformation parameters each time (e.g., translation L, P; translation L, S; rotation LR, AP; rotation LR, IS).

I will try this out.

You need to try different metrics, metric parameters, sampling percentage, maybe image proprocessing, etc. - until you see metric surface plots that have only one peak, at the optimal transformation parameters.

I will experiment with different metric parameters, sampling percentage etc.

Note that most modern registration frameworks, such as Elastix or ANTs come with default presets that should work well for most cases. So, if you don’t have special, particularly challenging images then they should work without the need for tuning any parameters.

I have been working with simpleITK but I will certainly try out Elastix for comparison purpose.

Thank you again.

8

Thank you for the insights.

The phantom contains a few major anatomical structures but they are made up of homogeneous materials. This kind of phantom is perfect for testing software workflows, visualization, application of transforms, etc. However, these images do not resemble real anatomical images in term of image content, details, and quality. You cannot use these data sets for evaluating automatic intensity-based image registration algorithms

That makes sense.
My aim is to design algorithm which should comply with " Report of the AAPM Radiation Therapy Committee Task Group No. 132". It should work with their dataset and results should be within tolerance specified in the report.

I will try out your suggestions and tune my algorithm further.
Thank you once again.

9

While the phantom data set cannot be used for assessing quality of the automatic intensity-based image registration algorithm, it can be used for testing most other aspects of the registration workflow.

1 Like