Registration Ct and MRI with multiresolution

drcjaudet · January 6, 2020, 3:16pm

Hello itk community,
My goals is to create a code in simple itk to perform quality contrôle of geometri distortion on MRI.
The key for following the evolution of the system performances is to register every new MRI with CT. The CT will be the standard as there is no geometric distortion of it. The image is a mahphan phantom.
CT:0.5x0.5x0.5mm^3
MRI:1x1x1mm^3

I did not succed to accuratly (less than 1mm) register the two image with the following code.
Can you give me some advice.

def multires_registration(fixed_image, moving_image, initial_transform, ImageSamplingPercentage):
registration_method = sitk.ImageRegistrationMethod()
registration_method.SetMetricAsMattesMutualInformation(numberOfHistogramBins=50)
registration_method.SetMetricSamplingStrategy(registration_method.RANDOM)
registration_method.SetMetricSamplingPercentage(float(ImageSamplingPercentage)/100)
registration_method.SetInterpolator(sitk.sitkLinear)
registration_method.SetOptimizerAsGradientDescent(learningRate=1.0, numberOfIterations=100, estimateLearningRate=registration_method.Once)
registration_method.SetOptimizerScalesFromPhysicalShift()
registration_method.SetInitialTransform(initial_transform)
registration_method.SetShrinkFactorsPerLevel(shrinkFactors = [4,2,1])
registration_method.SetSmoothingSigmasPerLevel(smoothingSigmas = [2,1,0])
registration_method.SmoothingSigmasAreSpecifiedInPhysicalUnitsOn()

final_transform = registration_method.Execute(fixed_image, moving_image)
print('Final metric value: {0}'.format(registration_method.GetMetricValue()))
print('Optimizer\'s stopping condition, {0}'.format(registration_method.GetOptimizerStopConditionDescription()))
return final_transform

and in the main code:

image_CT = sitk.Cast(image_CT, sitk.sitkFloat64)
image_IRM = sitk.Cast(image_IRM, sitk.sitkFloat64)
initial_transform = sitk.CenteredTransformInitializer(image_CT, image_IRM, sitk.Euler3DTransform(), sitk.CenteredTransformInitializerFilter.MOMENTS)
final_transform =multires_registration(image_CT, image_IRM, initial_transform, 5)
print("main: recalage  ok")

zivy · January 6, 2020, 5:59pm

Hello @drcjaudet,

Generally speaking the code looks reasonable, something went wrong in the copy paste of estimateLearningRate=registration_method.Once). I am not familiar with the “mahphan phantom” which I suspect is important in understanding the problem (assumptions/constraints imposed by the device).

To better understand the problem we need more details, beyond that you are not getting the desired accuracy:

Is the optimization not converging?
How are you judging the registration error:
a. Qualitatively, visual inspection with errors looking like they are greater than 1 voxel?
b. Quantitatively, you have a set of corresponding points localized in both images with sub-voxel accuracy (the diagonal of your MR voxel size, sqrt(3), is already close to twice the accuracy you want so the localization has to be sub-voxel).

If you are expecting a deformation (geometric distortion) in the MR, then using a rigid registration will obviously result in “errors” due to these deformations. If there is a sub-region of the MR image that does not have geometric distortions then you can use that region for rigid registration and the residual errors outside that region are the geometric distortion. Otherwise it is hard to separate between errors arising from the registration and those due to the geometric distortion.

drcjaudet · January 7, 2020, 8:41am

Dear Zivy,
thank’s for your answer. here some image of the phantom on CT and MRI. On the register image Ct is red and MRI green.
Ct

MRI
MRI

Register

the optimization is converging:
Final metric value: -0.246166961916
Optimizer’s stopping condition, GradientDescentOptimizerv4Template: Convergence checker passed at iteration 9.

a) I can see a rotation and translation (z) inaccuracy by eyes like on the picture.
b)

The program automatically compare the localisation of the center of the 158 small sphere in the Ct an MRI image. The difference of localisation directly measure the geometric distortion between the two image. It have to be less then 1mm.
Indeed i don’t expect geometrix distortion in the center of the image. I can use a mask on Ct but it will be hard to automatically create it on the new mri image, i guess?

Which estimateLearningRate will you use?
if you want i can provide you with the two images?

Have a nice day,
Cyril

zivy · January 7, 2020, 1:57pm

Hello @drcjaudet,

I don’t think that the issue is with the learning rate as the registration converges to what looks like a reasonable solution. Possibly what we are seeing are the effects of geometric distortion?

Below is a possible registration workflow that is based on the details you provided.

Assumptions
a. Central portion of the physical phantom exhibits minimal or no geometric distortion (it is close to the center of the image).
b. There is no error in the MR DICOM tags indicating patient orientation (see this Jupyter notebook for details).

Workflow

CenteredTransformInitializer.
Rigid registration using whole image.
Rigid registration using cropped CT, around the center of the phantom.

Each transformation that is output from the step goes as initial transform into the next step. The cropped CT is hard coded, and we are assuming that after step 2 the center of the phantom in both images is closely aligned (it appears to be). Step 3 is thus not influenced by the geometric distortions that are farther away from the region that isn’t distorted.

Finally, you should establish what is the best possible result you can get. I would localize the corresponding sphere centers in both CT and MR and use the LandmarkBasedTransformInitializer. This gives you the best result you can hope to achieve. If this isn’t accurate enough then chances are that automatic intensity based registration won’t do better.

drcjaudet · January 8, 2020, 3:19pm

Hello Zivy,

for step 2 it seems to work better with:
estimateLearningRate=registration_method.EachIteration

thank for the suggestion. I will try step 3 and let you know.
Thank you,
Cyril

drcjaudet · January 10, 2020, 1:31pm

Hello ITK,

i try the 3 step registration. The result seems a little bit more accurate but more robust. Even is the region of interest is not the same on the two modalities it work well.
In step 3 only the CT image is cropped. Do you think it will have an added value to also crop the mri on the same region?

If somenone is interested here is the code:

def multires_registration(fixed_image, moving_image, initial_transform, ImageSamplingPercentage):
    registration_method = sitk.ImageRegistrationMethod()
    registration_method.SetMetricAsMattesMutualInformation(numberOfHistogramBins=50)
    registration_method.SetMetricSamplingStrategy(registration_method.RANDOM)
    registration_method.SetMetricSamplingPercentage(float(ImageSamplingPercentage)/100)
    registration_method.SetInterpolator(sitk.sitkLinear)
    registration_method.SetOptimizerAsGradientDescent(learningRate=1.0, numberOfIterations=100, estimateLearningRate=registration_method.EachIteration) #Once
    registration_method.SetOptimizerScalesFromPhysicalShift() 
    registration_method.SetInitialTransform(initial_transform)
    registration_method.SetShrinkFactorsPerLevel(shrinkFactors = [4,2,1])
    registration_method.SetSmoothingSigmasPerLevel(smoothingSigmas = [2,1,0])
    registration_method.SmoothingSigmasAreSpecifiedInPhysicalUnitsOn()

    transform = registration_method.Execute(fixed_image, moving_image)
    print('Final metric value: {0}'.format(registration_method.GetMetricValue()))
    print('Optimizer\'s stopping condition, {0}'.format(registration_method.GetOptimizerStopConditionDescription()))
    return transform

and in main:

image_CT = sitk.Cast(image_CT, sitk.sitkFloat64)
    image_IRM = sitk.Cast(image_IRM, sitk.sitkFloat64)
    ####################recalage primaire de la sommes des moments des images###########
    initial_transform = sitk.CenteredTransformInitializer(image_CT, image_IRM, sitk.Euler3DTransform(), sitk.CenteredTransformInitializerFilter.MOMENTS)
    #####################recalage multiresolution rigide des images #################################
    medium_transform =multires_registration(image_CT, image_IRM, initial_transform, 5)
    #print("main: recalage  ok")
    #############################recalage multiresolution du centre du phantome moins soumis aux distortions geom#######################
    ####crop sur les 5 plus grosse sphéres#######################
    mask_label=sitk.BinaryThreshold(label_CT, 250, 254, 1, 0) #creer un masque sur les 5 plus grosse spheres
    stats= sitk.LabelIntensityStatisticsImageFilter()
    stats.Execute(mask_label, image_CT)
    delta=10 #additionnal voxel to crop volume to avoid border problem
    LowerBondingBox=[stats.GetBoundingBox(1)[0]-delta,stats.GetBoundingBox(1)[1]-delta,stats.GetBoundingBox(1)[2]-delta]
    UpperBondingBox=[image_CT.GetSize()[0]-(stats.GetBoundingBox(1)[0]+stats.GetBoundingBox(1)[3]+delta),image_CT.GetSize()[1]-(stats.GetBoundingBox(1)[1]+stats.GetBoundingBox(1)[4]+delta),image_CT.GetSize()[2]-(stats.GetBoundingBox(1)[2]+stats.GetBoundingBox(1)[5]+delta)]   
    image_CT_crop=cropImagefctLabel(image_CT, LowerBondingBox, UpperBondingBox ) 
    final_transform =multires_registration(image_CT_crop, image_IRM, medium_transform, 20)
    ##################reechantillonage############
    image_IRM_reg=reechantillonage(image_CT,image_IRM, final_transform)     
    su.PushToSlicer(image_IRM_reg, "verif_image_IRMrecaler",1) 
    #print time.ctime()
    print("main: reechantillonage  ok")

zivy · January 10, 2020, 3:10pm

Hello @drcjaudet,

Very nice of you to post the code for the benefit of other people that my be dealing with similar problems!

With respect to also cropping the MRI, I would advise against that, as you may hurt the registration if some portion of the CT is mapped outside the cropped MRI region while it is inside the original MRI volume. That missing information is relevant for registration. This is a theoretical answer, but without any real familiarity with your data.

When in doubt, try it out! If you try it out and find that cropping both volumes gives you more robust and accurate results you have empirical results from the world. All that remains is to identify the theoretical reason underlying your observation, what makes this specific data behave in such a way that generic data isn’t expected to do.

lassoan · January 14, 2020, 10:38pm

I have experienced many times that registration fails (converges to incorrect position) if one region covers a different region compared to the other but succeeds if the images are cropped to the common region (even if their initial alignment was good, regardless of which one was the fixed image). I suspect that it might be do to insufficient sampling but I did not investigate it much.