CT and PET coregistration with alternating results

Keyn34 · February 10, 2023, 10:36am

Hello everyone,

I am a bit confused about the registration results I got from performing a rigid alignment between a CT and a PET (of the same subject, which makes it even more confusing for me). I can run it once and get a proper alignment, but it might be off on the next run.

I think the best way to find the problem is to share the code:

def rigid(fixed_image, moving_image, aligned_image=None):
    # Define the registration method
    registration_method = SimpleITK.ImageRegistrationMethod()
    registration_method.SetInterpolator(SimpleITK.sitkLinear)

    # Inital alignment of the images
    initial_transform = SimpleITK.CenteredTransformInitializer(fixed_image, moving_image,
                                                               SimpleITK.Euler3DTransform(),
                                                               SimpleITK.CenteredTransformInitializerFilter.GEOMETRY)
    registration_method.SetInitialTransform(initial_transform, inPlace=True)

    # Only to check the inital alignment
    if aligned_image is not None:
        aligned_image_parent_directory = os.path.dirname(aligned_image)
        centered_image_path = os.path.join(aligned_image_parent_directory, 'rigid_centered_image.nii.gz')
        centered_resampled_image = SimpleITK.Resample(moving_image, fixed_image, initial_transform, SimpleITK.sitkLinear)
        SimpleITK.WriteImage(centered_resampled_image, centered_image_path)

    # Similarity metric settings
    registration_method.SetMetricAsMattesMutualInformation(numberOfHistogramBins=50)
    registration_method.SetMetricSamplingStrategy(registration_method.EachIteration)
    registration_method.SetMetricSamplingPercentage(0.01)

    # Optimizer settings
    registration_method.SetOptimizerAsGradientDescent(learningRate=2.0,
                                                      numberOfIterations=200,
                                                      convergenceMinimumValue=1e-6,
                                                      convergenceWindowSize=10)
    registration_method.SetOptimizerScalesFromPhysicalShift()

    # Multi-resolution framework settings
    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)

    # Always check the reason optimization terminated.
    print('Final metric value: {0}'.format(registration_method.GetMetricValue()))
    print('Optimizer\'s stopping condition, {0}'.format(registration_method.GetOptimizerStopConditionDescription()))

    if aligned_image is not None:
        resampled_image = SimpleITK.Resample(moving_image, fixed_image, final_transform)
        SimpleITK.WriteImage(resampled_image, aligned_image)

    return final_transform

Sometimes, the registration stops around iteration 9 and around iteration 37.
The initial alignment looks alignment after the initial transform always looks good, so that works fine.

I read through the registration notebooks (6* series) and tried to replicate everything as closely as possible, but obviously, there is something that I am missing and/or not understanding. Maybe you can hint me towards a specific section of the notebooks.

I don’t understand why I receive alternating results. I also tried it with affine registration, but there the optimizer reached the maximum number of iterations (200). So most likely, the problem might the optimizer that I set and the parameters for it. I also tried it with registration_method.Random first because I thought that the error might stem from there.

Any apparent mistakes in the code or the approach? I really appreciate any help you can provide.

zivy · February 10, 2023, 1:58pm

Hello @Keyn34,

First observation, wrong constant value given to the SetMetricSamplingStrategy:

registration_method.SetMetricSamplingStrategy(registration_method.EachIteration)

It expects one of the values from the MetricSamplingStrategyType which are NONE, REGULAR, RANDOM.

Not sure if it is lucky or unlucky, but we have:

registration_method.EachIteration == registration_method.RANDOM

So the registration is using a random sampling with a sample size of 1% of the image pixels (SetMetricSamplingPercentage(0.01)). This may work, but is a small sample size which can result in instability.

Some suggestions:

Change the SetMetricSamplingStrategy to RANDOM so using the expected strategy.
Comment out the multi-resolution framework settings (three lines) so that it is easy to debug. Possibly you may not need to go multi-resolution, always try that first before adding this complexity.
Most importantly, increase the sample size so that you get stable results. Try increasing to 0.1 and see if that increases stability. The specific sample size depends on your expected inputs, the required stability and the time constraints (more samples longer time).

Keyn34 · February 10, 2023, 10:13pm

Hello @zivy,

First, thank you for the explanation of the parameters. That also explains why I didn’t notice any difference between registration_method.EachIteration and registration_method.RANDOM.

I set the SetMetricSamplingStrategy to RANDOM now and set the sample size to 0.1. This resulted in stable results and reproducible registrations.

Something that confuses me, though: I also removed the multi-resolution framework, which sped up the rigid and affine registration I tried. It even made the affine registration perform way faster and more stable. It stops now after iteration 17 instead of 200, so it basically made it work. I always had the impression and understanding that a multi-resolution approach always speeds up the process and makes the registration more robust. Is there an explanation or a resource that might enlighten me on that?

Thank you so much again!

mihail.isakov · February 11, 2023, 1:39pm

Are you sure that fixed-moving approach is required? CT / PET should be aligned in physical space. You could check Frame of Reference UIDs of the series (if they were not messed by a bad anonymizer).
S. this @zivy post

The essential line is

# Resample PET onto CT grid using default interpolator and identity transformation.
pet_image_resampled = sitk.Resample(pet_image, ct_image)

Shalom_Rochman · May 1, 2024, 1:26pm

In certain scenarios, PET and CT scans may not align. For instance, in cases of liver metastatic lesions, a CT scan is initially performed without IV contrast. This CT scan serves as the basis for attenuation correction in the subsequent PET scan, which follows immediately. Following this procedure, the patient receives IV contrast, and approximately 10 minutes after the PET scan, a venous-phase CT scan is conducted. While the Frame of Reference UID may match between the venous-phase CT scan and the PET scan, the patient’s position has already shifted significantly. In such cases, a resampling strategy proves inadequate; image registration becomes necessary.