I figured out why I was getting a terrible landmark-aligned image:
When it comes to the use of LandmarkBasedTransformInitializer the numberOfControlPoints should be as low as possible (the BSpline order + 1), or 4 in my case.
Setting numberOfControlPoints to even moderately higher values highly distorts the image, and presummably that really bad starting point will lead the optimization astray.
I suppose this is quite obvious but I didn’t conscientiously set a high numberOfControlPoints for LandmarkBasedTransformInitializer. But since I wrapped it up in the registration method function, the same numberOfControlPoints that I wanted for the registration ended up being used for LandmarkBasedTransformInitializer. Hopefully others won’t make the same mistake.
The bad news is that I’m still getting the ComputeJacobianWithRespectToPosition error.
I’ve taken a slightly different approach this time. After succeeding in getting an affine registration with LandmarkBasedTransformInitializer and composing the LandmarkTx and OptimizedTx (example) in the ITKv4 framework (example), I made changes to it for a bspline transform (example):
def BsplineRegWithFiducials(FixIm, MovIm, FixPts, MovPts,
NumControlPts=8, SamplingPercentage=5,
LearningRate=5.0, NumIters=100):
import SimpleITK as sitk
import numpy as np
from ipywidgets import interact, fixed
dim = FixIm.GetDimension()
# Get landmark transform:
LandmarkTx = sitk.LandmarkBasedTransformInitializer(transform=sitk.BSplineTransform(dim),
fixedLandmarks=FixPts,
movingLandmarks=MovPts,
referenceImage=FixIm,
numberOfControlPoints=4)
SamplingPercentage = SamplingPercentage/100
ScaleFactors = [1, 2, 4]
MeshSize = [NumControlPts // ScaleFactors[-1]] * dim
OptimizedTx = sitk.BSplineTransformInitializer(image1=FixIm,
transformDomainMeshSize=MeshSize)
RegMethod = sitk.ImageRegistrationMethod()
#RegMethod.SetInitialTransform(InitialTx, True)
# Set the initial moving and optimized transforms in ITKv4
RegMethod.SetMovingInitialTransform(LandmarkTx)
RegMethod.SetInitialTransformAsBSpline(OptimizedTx, inPlace=True,
scaleFactors=ScaleFactors)
RegMethod.SetMetricAsMattesMutualInformation(numberOfHistogramBins=50)
RegMethod.SetMetricSamplingStrategy(RegMethod.RANDOM)
RegMethod.SetMetricSamplingPercentage(SamplingPercentage)
RegMethod.SetOptimizerAsGradientDescentLineSearch(learningRate=LearningRate,
numberOfIterations=NumIters,
convergenceMinimumValue=1e-4,
convergenceWindowSize=5)
RegMethod.SetOptimizerScalesFromPhysicalShift()
RegMethod.SetInterpolator(sitk.sitkLinear)
RegMethod.SetShrinkFactorsPerLevel([4, 2, 1])
RegMethod.SetSmoothingSigmasPerLevel([4, 2, 1])
RegMethod.AddCommand(sitk.sitkStartEvent, StartPlot)
RegMethod.AddCommand(sitk.sitkEndEvent, EndPlot)
RegMethod.AddCommand(sitk.sitkMultiResolutionIterationEvent, UpdateMultiresIterations)
RegMethod.AddCommand(sitk.sitkIterationEvent, lambda: PlotValues(RegMethod))
OptimizedTx = RegMethod.Execute(FixIm, MovIm)
# Compose the transformations:
FinalTx = sitk.CompositeTransform(OptimizedTx)
FinalTx.AddTransform(LandmarkTx)
#FinalTx = sitk.BSplineTransform(FinalTx)
print(f'Final metric value: {RegMethod.GetMetricValue()}')
print("Optimizer stop condition:",
f"{RegMethod.GetOptimizerStopConditionDescription()}")
# Resample MovIm using FinalTx to get the registered image:
RegIm = sitk.Resample(MovIm, FixIm, FinalTx, sitk.sitkLinear, 0.0,
MovIm.GetPixelID())
return RegIm
When I use the above method I get the following error:
Exception thrown in SimpleITK ImageRegistrationMethod_Execute: d:\a\1\sitk-build\itk-prefix\include\itk-5.1\itkImageToImageMetricv4GetValueAndDerivativeThreaderBase.hxx:268:
Exception in GetValueAndDerivativeProcessPoint:
d:\a\1\sitk-build\itk-prefix\include\itk-5.1\itkBSplineBaseTransform.h:302:
itk::ERROR: itk::ERROR: BSplineTransform(0000027AA7FCA850): ComputeJacobianWithRespectToPosition not yet implemented for BSplineTransform
which occurs at the line:
OptimizedTx = RegMethod.Execute(FixIm, MovIm)
Searching for any info on ComputeJacobianWithRespectToPosition has not led to any clues.
Question 1 (error):
Could someone please sanity check my method / offer any clues on the error / suggestions on how to resolve it?
Question 2 (no. of control points):
It’s not clear to me whether a different numberOfControlPoints for LandmarkTx and OptimizedTx is problemmatic. Should my selection of scaleFactors bridge the numberOfControlPoints in those two transforms?
I’ve not been able to find any examples where a composed transform was used that involved bsplines.
Question 3 (the transform used in LandmarkBasedTransformInitializer):
When I encountered the error I tried using sitk.VersorRigid3DTransform() or sitk.AffineTransform(dim) instead of sitk.BSplineTransform(dim) when defining LandmarkBasedTransformInitializer.
Doing so avoided the error, the registration proceeded, the plot of metric values v iterations looked “smooth” but when I plotted RegIm it was identical to AlignedIm.
I can confirm that the composite transform had a VersorRigid3Transform or AffineTransform after the BSplineTransform.
Can anyone comment on the merits of using a global v local transform for initializing a registration using a local transform? And more importantly, any ideas as to why there were no local deformations?
Question 4 (composing transforms):
Although I didn’t need to do it in the above method, I will later need to flatten the composition, which I thought was done with:
ComposedTx.FlattenTransform()
but when I ran that with the composed transform FinalTx used for the (successful) affine registration method using fiducials I got None.
Any ideas why?