BSpline registration for non-overlapping images

Algorithms
, , ,
1

Having run into problems trying to use landmarks to initiate a bspline registration of 3D CT and MR images, I decided to try a workaround: I cropped the CT image such that the field of view (in all three directions) is closely matched to that of the MR image stack.

And I also decided to make the registration a bit simpler as a sanity check. So I performed an affine (rather than BSpline) registration of the CT to MR image:

CTIF = sitk.CenteredTransformInitializerFilter.GEOMETRY
#CTIF = sitk.CenteredTransformInitializerFilter.MOMENTS

transf = sitk.AffineTransform(FixIm.GetDimension())

InitialTx = sitk.CenteredTransformInitializer(FixIm, MovIm, transf, CTIF)

RegMethod = sitk.ImageRegistrationMethod()
RegMethod.SetInitialTransform(InitialTx, inPlace=False)
RegMethod.SetMetricAsMattesMutualInformation(numberOfHistogramBins=50)
RegMethod.SetMetricSamplingStrategy(RegMethod.RANDOM)
RegMethod.SetMetricSamplingPercentage(0.05)

RegMethod.SetShrinkFactorsPerLevel(shrinkFactors=[4,2,1])
RegMethod.SetSmoothingSigmasPerLevel(smoothingSigmas=[2,1,1])
RegMethod.SmoothingSigmasAreSpecifiedInPhysicalUnitsOn()

RegMethod.SetInterpolator(sitk.sitkLinear)

RegMethod.SetOptimizerAsGradientDescent(learningRate=1.0, 
                                        numberOfIterations=500, 
                                        estimateLearningRate=RegMethod.Once)
RegMethod.SetOptimizerScalesFromPhysicalShift()

FinalTx = RegMethod.Execute(FixIm, MovIm)

RegIm = sitk.Resample(MovIm, FixIm, FinalTx, sitk.sitkLinear, 0.0, 
                      MovIm.GetPixelID())

using either sitk.CenteredTransformInitializerFilter.GEOMETRY or sitk.CenteredTransformInitializerFilter.MOMENTS, I get better or worse results, but they both seem reasonable.

However I’m really interested in performing a deformable registration (following this example, and these three bspline examples: example1, example2, example3), but unlike sitk.CenteredTransformInitializer, sitk.BSplineTransformInitializer doesn’t appear to allow for use of GEOMETRY or MOMENTS to help guide the initial alignment.

Despite having cropped the CT image to roughly match the field-of-view of the MR image, the z-coordinates of the two images are offset by around 750 mm. So of course if I attempt to run the BSpline registration (see code further below) I get:

itk::ERROR: MattesMutualInformationImageToImageMetricv4: Joint PDF summed to zero

because there is no overlap between the two images.

I’ve seen it written multiple times that when it comes to initializations there are generally four options:

  1. Do nothing (–> “Joint PDF summed to zero” error)
  2. CenteredTransformInitializer using GEOMETRY or MOMENTS (–> not an option for BSplineTransformInitializer)
  3. Exhaustive optimizer (–> still use of CenteredTransformInitializer)
  4. LandmarkBasedTransformInitializer (–> “MultiThreader: Exception occurred during SingleMethodExecute” error)

I feel that there must be something I’ve missed out on, or have misunderstood in applying BSpline registrations. Any insight would be much appreciated.

Following is the code used to attempt the bspline registration:

GridSpacing = [50.0, 50.0, 50.0]
ImPhySize = [size*spacing for size, spacing in zip(FixIm.GetSize(), 
                                                   FixIm.GetSpacing())]
MeshSize = [int(sz/sp + 0.5) for sz, sp in zip(ImPhySize, GridSpacing)]
MeshSize = [int(sz/4 + 0.5) for sz in MeshSize]

InitialTx = sitk.BSplineTransformInitializer(image1=FixIm, 
                                             transformDomainMeshSize=MeshSize, 
                                             order=3)

RegMethod = sitk.ImageRegistrationMethod()
RegMethod.SetInitialTransformAsBSpline(InitialTx, inPlace=False,
                                       scaleFactors=[1,2,4])
RegMethod.SetMetricAsMattesMutualInformation(numberOfHistogramBins=50)
RegMethod.SetMetricSamplingStrategy(RegMethod.RANDOM)
RegMethod.SetMetricSamplingPercentage(0.05)

RegMethod.SetShrinkFactorsPerLevel(shrinkFactors=[4,2,1])
RegMethod.SetSmoothingSigmasPerLevel(smoothingSigmas=[2,1,1])
RegMethod.SmoothingSigmasAreSpecifiedInPhysicalUnitsOn()

RegMethod.SetInterpolator(sitk.sitkLinear)

RegMethod.SetOptimizerAsLBFGS2(solutionAccuracy=1e-2, 
                               numberOfIterations=100, 
                               deltaConvergenceTolerance=0.01)
RegMethod.SetOptimizerScalesFromPhysicalShift()

FinalTx = RegMethod.Execute(FixIm, MovIm)

RegIm = sitk.Resample(MovIm, FixIm, FinalTx, sitk.sitkLinear, 0.0, 
                      MovIm.GetPixelID())