Registering vertebra CT and MRI

[Ray_Jian]

Hello, All. I am trying to register vertebra in spine CT and MRI.

The average size of the vertebra MRI is (200-300, 100-200, 15).
The average number of intensity in MRI is around 500. Which is calculated by this.

min_max = np.percentile(new_scan_array.flatten(),[2,98])
num_bins = int(min_max[1]-min_max[0])

The average size of the vertebra CT is (200-300, 100-200, 30-50).
THe average number of intensity in CT is around 1000.

vertebra 19:

vertebra 9:

Here are two examples of how the crop MRI and CT looks like. (The examples shown are T2-weighted, I have also tried with T1-weighted)

For Initialization, I have used the centorid of the interested vertebra. (The centroids are marked manually in the whole spine scan, and when I cropped the CT and MRI, I maintained the relative position of all vertebrae by resetting the origin.)

initial_transform = sitk.LandmarkBasedTransformInitializer(sitk.VersorRigid3DTransform(), 
                                                            fixed_image_points_flat, 
                                                            moving_image_points_flat)

Here is the visualization of the initialization:
vertebra 19:

vertebra 9:

Then I set up the registration method:

registration_method = sitk.ImageRegistrationMethod()

# Similarity metric settings.
registration_method.SetMetricAsMattesMutualInformation(numberOfHistogramBins=1000)
# registration_method.SetMetricAsANTSNeighborhoodCorrelation(radius=20)
registration_method.SetMetricSamplingStrategy(registration_method.NONE)
# registration_method.SetMetricSamplingPercentage(0.1)

registration_method.SetInterpolator(sitk.sitkLinear)

# Optimizer settings.
registration_method.SetOptimizerAsGradientDescent(learningRate=1, numberOfIterations=100, 
                                                  convergenceMinimumValue=1e-6, convergenceWindowSize=10)
# registration_method.SetOptimizerAsRegularStepGradientDescent(learningRate=1,minStep = 1e-6, relaxationFactor = 0.9, 
#                                                                  numberOfIterations=100)

registration_method.SetOptimizerScalesFromPhysicalShift()

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

I used MattesMutualInformation and GradientDescent, I set the numberOfHistogramBins to 1000 because the number of intensity values is around 1000. (Please correct me if my intuition is wrong here.)
The numberOfIteration I tired is 100 or 200.

Here are the curves of the metric value:
Vertebra 19:

Vertebra 9:

And here is the final registration result:
Vertebra 19:

Vertebra 9:

As you can see, vertebra 9 is registered well while vertebra 19 is not so good. I want the algorithm to be robust to every vertebra, so do you have any ideas for the method or parameters that I have used? It would be very nice of you to provide some practical advice. If there is more needed to be seen, please inform me.

By the way, I have also tried with the ITK_reigstration_v4 method.

optimized_transform = sitk.VersorRigid3DTransform()    
registration_method.SetMovingInitialTransform(initial_transform)
registration_method.SetInitialTransform(optimized_transform, inPlace=False)
final_transform_v4 = registration_method.Execute(fixed_image, moving_image)
final_transform_v4.AddTransform(initial_transform)

But it didn’t work at all, the result made no sense.

Best.

[dzenanz (Dženan Zukić)]

@Ray_Jian

numberOfHistogramBins to 1000

This probably defeats the method. It should be a much smaller number. Try not setting it at all, so default value is used. CT-MRI is most common multi-modality registration pair.

[Ray_Jian]

Thanks a lot!! That’s working! Is there any advice on other parameters?

[Ray_Jian]

By the way, if I want to design my own similarity metric in ITK (cpp) and use it in SimpleITK (python), what should I do? Is there any tutorial for that?
I want to use LC2 metric (http://campar.in.tum.de/Main/LC2Code),
but it is not implemented in SimpleITK, so I was thinking that I may could add this in SimpleITK and use the GradientDescent optimizer to do registration.
Thanks a lot for your kind answer!

[zivy (Ziv Yaniv)]

Hello @Ray_Jian,

Unfortunately, we don’t have a tutorial on how to add your own metric. We’re always looking for contributions like that (would go on the toolkits read-the-docs site).

To add a metric to SimpleITK you first need to implement it in ITK. It is recommended that you start with the closest existing metric (see Metricsv4).

You then need to add it to SimpleITK, edit the sitkImageRegistrationMethod_CreateMetric.hxx file and the sitkImageRegistrationMethod.h file.

As this is a custom SimpleITK version you will need to build it from source.

[dzenanz (Dženan Zukić)]

Also take a look at a short answer and a longer discussion.

[lassoan (Andras Lasso)]

Medical image registration using classic (not machine learning based) methods is a very mature field. Thousands of researchers spent decades developing registration algorithms, and most solvable problems got solved. So, trying to come up with a new metric that beats the existing ones would be hard. Trying existing methods, tune and maybe combine them would be probably a better investment of your time.

For example, you can try the ITK-based Elastix (or SimpleElastix) toolkit. It adds many new metrics and optimization strategies.

If you must develop new algorithms (e.g., because you need to demonstrate that you can develop new computational algorithms to graduate) then you may want to study machine learning based methods. Not because ML-based methods are better, but because they are new, there is much more to explore, and you would obtain skills that are highly sought after nowadays.