Query on Regularization with BSplineTransform and GradientDescent in SimpleITK

Beginner Questions
, ,
1

Hello everyone,

I am currently working with SimpleITK for image registration, specifically employing sitk.BSplineTransform and the SetOptimizerAsGradientDescent optimizer.

I’ve been reviewing a Python example that demonstrates B-spline registration (provided below for reference).

My observation is that, in this code, there aren’t any explicit options or parameters for incorporating a regularization term directly within the SetOptimizerAsGradientDescentLineSearch method or when configuring the registration process. While SetShrinkFactorsPerLevel and SetSmoothingSigmasPerLevel can provide a form of implicit smoothing across multi-resolution levels, I’m looking for a more explicit regularization mechanism.

My question is: Is it possible to directly incorporate a regularization term (e.g., a penalty for excessive deformation or a smoothness prior for the B-spline coefficients) when using sitk.GradientDescent with sitk.BSplineTransform in SimpleITK? If so, could you please advise on how this can be implemented or configured?

Any insights, documentation references, or illustrative examples would be greatly appreciated.

Thank you

#!/usr/bin/env python

from __future__ import print_function

import SimpleITK as sitk
import sys
import os


def command_iteration(method) :
    print("{0:3} = {1:10.5f}".format(method.GetOptimizerIteration(),
                                     method.GetMetricValue()))
    print("\t#: ", len(method.GetOptimizerPosition()))


def command_multi_iteration(method) :
    print("--------- Resolution Changing ---------")


if len ( sys.argv ) < 4:
    print( "Usage: {0} <fixedImageFilter> <movingImageFile> <outputTransformFile>".format(sys.argv[0]))
    sys.exit ( 1 )


fixed = sitk.ReadImage(sys.argv[1], sitk.sitkFloat32)

moving = sitk.ReadImage(sys.argv[2], sitk.sitkFloat32)

transformDomainMeshSize=[10]*moving.GetDimension()
tx = sitk.BSplineTransformInitializer(fixed,
                                      transformDomainMeshSize )

print("Initial Parameters:");
print(tx.GetParameters())

R = sitk.ImageRegistrationMethod()
R.SetMetricAsMattesMutualInformation(50)
R.SetOptimizerAsGradientDescentLineSearch(5.0, 100,
                                          convergenceMinimumValue=1e-4,
                                          convergenceWindowSize=5)
R.SetOptimizerScalesFromPhysicalShift( )
R.SetInitialTransform(tx)
R.SetInterpolator(sitk.sitkLinear)

R.SetShrinkFactorsPerLevel([6,2,1])
R.SetSmoothingSigmasPerLevel([6,2,1])

R.AddCommand( sitk.sitkIterationEvent, lambda: command_iteration(R) )
R.AddCommand( sitk.sitkMultiResolutionIterationEvent, lambda: command_multi_iteration(R) )

outTx = R.Execute(fixed, moving)

print("-------")
print(outTx)
print("Optimizer stop condition: {0}".format(R.GetOptimizerStopConditionDescription()))
print(" Iteration: {0}".format(R.GetOptimizerIteration()))
print(" Metric value: {0}".format(R.GetMetricValue()))

sitk.WriteTransform(outTx,  sys.argv[3])

if ( not "SITK_NOSHOW" in os.environ ):

    resampler = sitk.ResampleImageFilter()
    resampler.SetReferenceImage(fixed);
    resampler.SetInterpolator(sitk.sitkLinear)
    resampler.SetDefaultPixelValue(100)
    resampler.SetTransform(outTx)

    out = resampler.Execute(moving)
    simg1 = sitk.Cast(sitk.RescaleIntensity(fixed), sitk.sitkUInt8)
    simg2 = sitk.Cast(sitk.RescaleIntensity(out), sitk.sitkUInt8)
    cimg = sitk.Compose(simg1, simg2, simg1//2.+simg2//2.)
    sitk.Show( cimg, "ImageRegistration1 Composition" )
2

Hello @blue_sky,

This is where you are starting to transition from SimpleITK to ITK or possibly ANTs. SimpleITK does not support explicit regularization when using the BSplineTransform. ITK does via the BSplineSmoothingOnUpdateDisplacementFieldTransform class. The ANTs folks contributed that code to ITK, see the “Explicit B-spline regularization in diffeomorphic image registration” paper. I believe the relevant Python documentation is here.

@ntustison, please correct me if I’m wrong and possibly provide some more insight :wink:

1 Like
3

The BSplineSSmoothingOnUpdateDisplacementFieldTransform and GaussianSmoothingOnUpdateDisplacementFieldTransform are available option for the SimpleITK DisplacementFieldTransform: SimpleITK: itk::simple::DisplacementFieldTransform Class Reference

1 Like
4

Hey all,

Currently away and I’ll be sure to provide a more complete answer when I get back but my guess is that @blue_sky would be more interested in what the elastix folk have done in terms of external regularization (e.g., Itk::Transforms supporting spatial derivatives). I don’t know how much has been ported to ITK or SimpleITK proper but I suspect it’s a good place to start exploring.

Nick

5

In case people are interested about what’s in ITK/SimpleITK, specifically the classes that @zivy and @blowekamp point to, this is completely separate work from the explicit regularization terms that were proposed in Dan Rueckert’s original FFD paper and later developed in Elastix.

In addition to the paper mentioned by @zivy, where we describe our development of the B-spline variant of Symmetric Normalization, I wrote an earlier paper, Directly manipulated free-form deformation image registration, that discusses the relation of classic B-spline image registration (without explicit regularization) to the Demons algorithm where one can simply substitute Gaussian smoothing in the latter with a specific type of B-spline smoothing which turns out to be equivalent to a type of preconditioner on the former. This is what is contained in ITK/SimpleITK.

In addition to the implicit smoothing of B-splines, SyN’s inverse constraint as well as the general preferred use of a localized neighborhood cross-correlation similarity metric also add regularization to the final displacement field(s). Coupling these considerations with the computational penalty associated with inclusion of additional regularization terms, we haven’t prioritized their development in ANTs.

2 Likes