How to initialize CenteredTransformInitializer for 3D-3D image registration for exhaustive search optimizer

Hi All,

I’m following the 2D image registration example(python) for exhaustive search optimizer: https://examples.itk.org/src/numerics/optimizers/exhaustiveoptimizer/documentation

I want to do the same for 3D-3D image registration. But I couldn’t intialize CenteredTransformInitializer with itk.VersorRigid3DTransform[itk.D]. Could you please share a sample code snippet.

This is my code:

    dimension = fixed_image.GetImageDimension()
     
    # Use cast filter to convert UC type images to SS
    InputImageType = itk.Image[itk.UC, dimension]
    OutputImageType = itk.Image[itk.SS, dimension]
    cast_filter = itk.CastImageFilter[InputImageType, OutputImageType].New()
    cast_filter.SetInput(fixed_image)
    cast_filter.Update()
    fixed_image = cast_filter.GetOutput()
    cast_filter.SetInput(moving_image)
    cast_filter.Update()
    moving_image = cast_filter.GetOutput()
    
    FixedImageType = type(fixed_image)
    MovingImageType = type(moving_image)

    if dimension == 3:
        transform = itk.Euler3DTransform[itk.D].New()
    else:
        raise Exception(f"Unsupported dimension: {dimension}")

    TransformInitializerType = itk.CenteredTransformInitializer[
        itk.VersorRigid3DTransform[itk.D],
        # itk.MatrixOffsetTransformBase[itk.D, dimension, dimension],
        FixedImageType,
        MovingImageType,
    ]
    initializer = TransformInitializerType.New(
        Transform=transform,
        FixedImage=fixed_image,
        MovingImage=moving_image,
    )
    initializer.InitializeTransform()

Thank you!

I was able to intialize CenteredTransformInitializer with VersorRigid3DTransform with following is the code snippet:

    dimension = fixed_image.GetImageDimension()
     
    # Use cast filter to convert UC type images to D
    InputImageType = itk.Image[itk.UC, dimension]
    OutputImageType = itk.Image[itk.D, dimension]
    cast_filter = itk.CastImageFilter[InputImageType, OutputImageType].New()
    cast_filter.SetInput(fixed_image)
    cast_filter.Update()
    fixed_image = cast_filter.GetOutput()
    cast_filter.SetInput(moving_image)
    cast_filter.Update()
    moving_image = cast_filter.GetOutput()
    
    FixedImageType = type(fixed_image)
    MovingImageType = type(moving_image)

    if dimension == 3:
        transform = itk.VersorRigid3DTransform[itk.D].New()
        # transform = itk.Euler3DTransform[itk.D].New()
    else:
        raise Exception(f"Unsupported dimension: {dimension}")

    TransformInitializerType = itk.CenteredTransformInitializer[
        itk.VersorRigid3DTransform[itk.D],
        # itk.Euler3DTransform[itk.D],
        # itk.MatrixOffsetTransformBase[itk.D, dimension, dimension],
        FixedImageType,
        MovingImageType,
    ]
    initializer = TransformInitializerType.New(
        Transform=transform,
        FixedImage=fixed_image,
        MovingImage=moving_image,
    )
    
    # initializer.MomentsOn()
    initializer.GeometryOn()
    initializer.InitializeTransform()
    
    if similarity_metric == "CC":
        metric = itk.CorrelationImageToImageMetricv4[
            FixedImageType, MovingImageType
        ].New()
    elif similarity_metric == "MI":
        metric = itk.MattesMutualInformationImageToImageMetricv4[
            FixedImageType, MovingImageType
        ].New()
    elif similarity_metric == "MSE":
        metric = itk.MeanSquaresImageToImageMetricv4[
            FixedImageType, MovingImageType
        ].New()


    optimizer = itk.ExhaustiveOptimizerv4[itk.D].New()

    def print_iteration():
        print(
            f"Iteration:"
            f"{optimizer.GetCurrentIteration()}\t"
            f"{list(optimizer.GetCurrentIndex())} \t"
            f"{optimizer.GetCurrentValue():10.4f}\t"
            f"{list(optimizer.GetCurrentPosition())}\t"
        )

    optimizer.AddObserver(itk.IterationEvent(), print_iteration)

    angles = 1
    optimizer.SetNumberOfSteps([int(angles / 2), 1, 1])
    # optimizer.SetStepLength(0.1)

    # Initialize scales and set back to optimizer
    scales = optimizer.GetScales()
    scales.SetSize(6)
    scales.SetElement(0, 1) #2.0 * np.pi / angles)
    scales.SetElement(1, 1) #2.0 * np.pi / angles)
    scales.SetElement(2, 1) #2.0 * np.pi / angles)
    scales.SetElement(3, 1.0)
    scales.SetElement(4, 1.0)
    scales.SetElement(5, 1.0)
    optimizer.SetScales(scales)  # Set the number of subdivisions per step length

    RegistrationType = itk.ImageRegistrationMethodv4[FixedImageType, MovingImageType]
    registration = RegistrationType.New(
        Metric=metric,
        Optimizer=optimizer,
        FixedImage=fixed_image,
        MovingImage=moving_image,
        InitialTransform=transform,
        NumberOfLevels=1,
    )

    try:
        registration.Update()

        print(
            f"MinimumMetricValue: {optimizer.GetMinimumMetricValue():.4f}\t"
            f"MaximumMetricValue: {optimizer.GetMaximumMetricValue():.4f}\n"
            f"MinimumMetricValuePosition: {list(optimizer.GetMinimumMetricValuePosition())}\t"
            f"MaximumMetricValuePosition: {list(optimizer.GetMaximumMetricValuePosition())}\n"
            f"StopConditionDescription: {optimizer.GetStopConditionDescription()}\t"
        )        

    except Exception as e:
        print(f"Exception caught: {e}")
        return None

But the optimizer values are high with CrossCorrelation metric:
Iteration:1220 [0.0, 2.0, 1.0, 36.0, 1.0, 0.0] 179769313486231570814527423731704356798070567525844996598917476803157260780028538760589558632766878171540458953514382464234321326889464182768467546703537516986049910576551282076245490090389328944075868508455133942304583236903222948165808559332123348274797826204144723168738177180919299881250404026184124858368.0000 [0.0, 0.9999999999, 0.0, -13.0, -179945855.0, -140706685389008.0]

This is a sample optimizer value for Mean Square metric:

Iteration:59    [0.0, 2.0, 1.0, 6.0, 0.0, 0.0]  179769313486231570814527423731704356798070567525844996598917476803157260780028538760589558632766878171540458953514382464234321326889464182768467546703537516986049910576551282076245490090389328944075868508455133942304583236903222948165808559332123348274797826204144723168738177180919299881250404026184124858368.0000  [0.0, 0.9999999999, 0.0, -939.0, -140194191693648.0, -140194191706992.0]
^CWARNING: In /work/ITK-source/ITK/Modules/Numerics/Optimizersv4/include/itkObjectToObjectMetric.hxx, line 575
MeanSquaresImageToImageMetricv4 (0x58990f0): No valid points were found during metric evaluation. For image metrics, verify that the images overlap appropriately. For instance, you can align the image centers by translation. For point-set metrics, verify that the fixed points, once transformed into the virtual domain space, actually lie within the virtual domain.

I’m getting error with Mattes Mutual Information metric as my images are not overlapped significantly:

Exception caught: /work/ITK-source/ITK/Modules/Registration/Metricsv4/include/itkMattesMutualInformationImageToImageMetricv4.hxx:312:
ITK ERROR: MattesMutualInformationImageToImageMetricv4(0xb5e9200): All samples map outside moving image buffer. The images do not sufficiently overlap. They need to be initialized to have more overlap before this metric will work. For instance, you can align the image centers by translation.

Could you please advice how to overcome this issues as I’m very new to ITK.

Thank you!

I have fixed the issue. I forgot to set number of steps for all 6 parameters when converting 2D code to 3D.

Thank you!

3 Likes