Using LandmarkBasedTransformInitializer for BSpline registrations

Algorithms
, ,
1

Hi all,

I originally posted this question on SimpleITK’s GitHub issues page, but became reminded about the Discourse forum while watching a SITK tutorial (somehow I completely forgot that the forum exists!). I’m re-posting it here hoping to reach a wider audience (and hopefully someone who’s done this before and knows where I’m going wrong).

I am aware of sitk.BSplineTransformInitializer(image1, transformDomainMeshSize, order) for performing BSpline registrations, but I am trying to initialize the registration using LandmarkBasedTransformInitializer.

Following is how I attempt to obtain a landmark-based transform for a ‘rigid’, ‘affine’ or ‘bspline’ transforms:

def GetLandmarkTx(FixIm, MovIm, FixPts, MovPts, 
                  Transform='rigid', NumControlPts=8):
    import SimpleITK as sitk
    
    if Transform == 'rigid':
        SitkTransform = sitk.VersorRigid3DTransform()
    elif Transform == 'affine':
        SitkTransform = sitk.AffineTransform(FixIm.GetDimension())
    elif Transform == 'bspline':
        SitkTransform = sitk.BSplineTransform(FixIm.GetDimension())
    else:
        msg = "Transform must be either 'rigid', 'affine' or 'bspline'."
        raise Exception(msg)
    
    if Transform == 'bspline':
        LandmarkTx = sitk.LandmarkBasedTransformInitializer(transform=SitkTransform, 
                                                            fixedLandmarks=FixPts, 
                                                            movingLandmarks=MovPts,
                                                            referenceImage=FixIm,
                                                            numberOfControlPoints=NumControlPts)
    else:
        LandmarkTx = sitk.LandmarkBasedTransformInitializer(transform=SitkTransform, 
                                                           fixedLandmarks=FixPts, 
                                                           movingLandmarks=MovPts,
                                                           referenceImage=FixIm)

    return LandmarkTx

The above works fine for ‘rigid’ and ‘affine’, but the following error occurs for ‘bspline’:

Exception thrown in SimpleITK LandmarkBasedTransformInitializer: C:...itkMultiThreader.cxx:395:
itk::ERROR: MultiThreader: Exception occurred during SingleMethodExecute

Unfortunately I’ve not been able to find any examples that use LandmarkBasedTransformInitializer for non-rigid registrations, so I’m not sure if I’m taking the correct approach. Any advice would be much appreciated.

1 Like
2

I stumbled across an example of the use of LandmarkBasedTransformInitializer for BSplineTransform and tried the following:

MeshSize = [NumControlPts] * MovIm.GetDimension()
    
InitTx = sitk.BSplineTransformInitializer(image1=FixIm,
                                          transformDomainMeshSize=MeshSize)
        
LandmarkTx = sitk.LandmarkBasedTransformInitializerFilter()
LandmarkTx.SetFixedLandmarks(FixPts)
LandmarkTx.SetMovingLandmarks(MovPts)
LandmarkTx.SetBSplineNumberOfControlPoints(NumControlPts)
LandmarkTx.SetReferenceImage(FixIm)
InitialTx = LandmarkTx.Execute(InitTx)

although I believe that is equivalent to:

dim = FixIm.GetDimension()
        
InitialTx = sitk.LandmarkBasedTransformInitializer(transform=sitk.BSplineTransform(dim),
                                                   fixedLandmarks=FixPts, 
                                                   movingLandmarks=MovPts,
                                                   referenceImage=FixIm,
                                                   numberOfControlPoints=NumControlPts)

I’m getting the same error using both methods, but the error is different than the one described in my original post (I have since updated SimpleITK from v1.2 to v2.0). The error is now:

itk::ERROR: itk::ERROR: BSplineScatteredDataPointSetToImageFilter(00000169F9BF0FC0): The reparameterized point component 5.07517 is outside the corresponding parametric domain of [0, 5).

I found this post about the same error, in which someone suggested to try a slightly larger spacing, but I’m not sure what was meant by that - a larger grid spacing?

I tried reducing the number of control points (to reduce the grid size and hence increase the grid spacing) but all that did was change the reparameterized point component and parametric domain - the point was still outside the domain.

I can’t see anything in the description of LandmarkBasedTransformInitializer that allows for changing the domain.

I’m out of ideas on how to overcome this dilemma. Does anyone have any suggestions?

3

Hello,

From your original code, the parameter FixImg defines the domain of the BSplineTransform transform. It looks like one of the points are outside this domain.

Brad

4

@blowekamp: Thanks for the suggestion.

I was able to perform an affine registration using LandmarkBasedTransformInitializer:

InitialTx = sitk.LandmarkBasedTransformInitializer(sitk.VersorRigid3DTransform(), 
                                                   FixPts, MovPts)

RegMethod = sitk.ImageRegistrationMethod()
RegMethod.SetInitialTransform(InitialTx, inPlace=False)
...

without receiving that error, and the result was good, i.e. MovIm (a full-body CT image) was successfully registered to FixIm (MR image of much limited extent) - a result that was not possible using CenteredTransformInitializer.

In the following plot are the fiducials that I used for FixIm (MR501) and MovIm (CT3).

20210426_131506_STS_004_fiducials
20210426_131506_STS_004_fiducials1871×8971 2.67 MB

Does that rule out an issue with my fiducials (such as there being at least one point in FixPts that lies beyond the extent of FixIm) or have I misunderstood your suggestion?

Thanks.

5

The affine based transforms are global. The BSplineTransform has a finite domain. So it does not rule it out.

7

@blowekamp: Thanks again.

I eliminated any pair of fiducials for which at least one lied on any boundary of its corresponding image (none were outside). I re-tried the registration but got the exact same error - in fact even the “reparameterized point component” was as before: 5.07517.

Then I excluded more aggressively - keeping only fiducials that were at least 4 pixels from the boundary (I found the limit to be 2 pix). That enabled the registration to proceed but the computational time was very long (74 min).

It turned out that TransformDomainMeshSize somehow became [20, 20, 20]. This was despite having set the number of control points to 8:

NumControlPts = 8

MeshSize = [NumControlPts] * MovIm.GetDimension()
    
InitTx = sitk.BSplineTransformInitializer(FixIm, MeshSize)

LandmarkTx = sitk.LandmarkBasedTransformInitializerFilter()
LandmarkTx.SetFixedLandmarks(FixPts)
LandmarkTx.SetMovingLandmarks(MovPts)
LandmarkTx.SetBSplineNumberOfControlPoints(NumControlPts)
LandmarkTx.SetReferenceImage(FixIm)
InitialTx = LandmarkTx.Execute(InitTx)

RegMethod = sitk.ImageRegistrationMethod()
RegMethod.SetInitialTransform(InitialTx, True)

where FixIm = MrIm and MovIm = CtIm.

(It’s also surprising given that by default BSplineNumberOfControlPoints = 4 for LandmarkBasedTransformInitializerFilter, and that’s what LandmarkTx was set to at the point of defining InitialTx.)

Is this expected behaviour or might it be a bug?

Anyhow, I added the following lines before RegMethod.SetInitialTransform(InitialTx, True):

InitialTx = sitk.BSplineTransform(InitialTx)
InitialTx.SetTransformDomainMeshSize(MeshSize)

which seemed to do the trick, and the registration took a much more reasonable 5 min.

Unfortunately the resulting registration was horrible. It seems that the problem lies with InitialTx. When I used it to resample MovIm I was expecting to get a result that looked as if CtIm had been cropped to a similar extent to that of MrIm. But the resulting image was very odd.

I repeated the resampling using the previous InitialTx (the one with [20, 20, 20] grid size) but it too looked very odd, albeit different from the other result.

I suspect that my approach of using LandmarkBasedTransformInitializerFilter is not necessarily what I want, since the setting of SetBSplineNumberOfControlPoints suggests that the transformation will be deformable (section 3.17 Point Set Registration in the sitk manual v5).

So I tried a different approach - I created LandmarkTx as I did when for affine registrations:

LandmarkTx = sitk.LandmarkBasedTransformInitializer(sitk.VersorRigid3DTransform(), 
                                                    FixPts, MovPts)

and

AlignedIm = sitk.Resample(MovIm, FixIm, AlignTx, sitk.sitkLinear, 0.0, 
                          MovIm.GetPixelID())

gave me sensible results (AlignedIm looked like CtIm was cropped to fit the extent of MrIm).

Then I created an initial BSpline Transform (as I do when performing a BSpline registration without fiducials):

MeshSize = [NumControlPts] * FixIm.GetDimension()

InitialTx = sitk.BSplineTransformInitializer(FixIm, MeshSize)

and combined both into a composite transform:

CompositeTx = sitk.Transform(InitialTx)
CompositeTx.AddTransform(LandmarkTx)

print(CompositeTx)

itk::simple::CompositeTransform
CompositeTransform (00000169F02F7920)
RTTI typeinfo: class itk::CompositeTransform<double,3>
Reference Count: 1
Modified Time: 13345167
Debug: Off
Object Name:
Observers:
none
Transforms in queue, from begin to end:

BSplineTransform (00000169FB1620A0)
RTTI typeinfo: class itk::BSplineTransform<double,3,3>
Reference Count: 1
Modified Time: 13345153
Debug: Off
Object Name:
Observers:
none
CoefficientImage: [ 00000169F2EF6A50, 00000169F2EF61E0, 00000169F2EF5F10 ]
TransformDomainOrigin: [-220.461, -137.269, -102.056]
TransformDomainPhysicalDimensions: [260.253, 260.253, 181.5]
TransformDomainDirection: 0.998543 -0.0277357 -0.0462969
0.0277059 0.999615 -0.0012846
0.0463148 0 0.998927

 TransformDomainMeshSize: [8, 8, 8]
 GridSize: [11, 11, 11]
 GridOrigin: [-250.992, -170.661, -126.226]
 GridSpacing: [32.5316, 32.5316, 22.6875]
 GridDirection: 0.998543 -0.0277357 -0.0462969

0.0277059 0.999615 -0.0012846
0.0463148 0 0.998927

VersorRigid3DTransform (00000169F648B180)
RTTI typeinfo: class itk::VersorRigid3DTransform
Reference Count: 1
Modified Time: 13345163
Debug: Off
Object Name:
Observers:
none
Matrix:
0.995695 0.0837221 0.0397639
-0.0847359 0.996101 0.0245329
-0.0375549 -0.0277968 0.998908
Offset: [4.33817, -7.13875, 730.617]
Center: [-82.026, 3.47923, -24.1324]
Translation: [4.02295, -0.793809, 733.627]
Inverse:
0.995695 -0.0847359 -0.0375549
0.0837221 0.996101 -0.0277968
0.0397639 0.0245329 0.998908
Singular: 0
Versor: [ -0.0130977, 0.0193522, -0.0421635, 0.998837 ]
End of MultiTransform.
<<<<<<<<<<
TransformsToOptimizeFlags, begin() to end():
0 1
TransformsToOptimize in queue, from begin to end:
End of TransformsToOptimizeQueue.
<<<<<<<<<<
End of CompositeTransform.
<<<<<<<<<<

Then carrying on with the registration method:

RegMethod = sitk.ImageRegistrationMethod()
RegMethod.SetInitialTransform(CompositeTx, True)

# setting of metrics, optimizer, etc. not shown

FinalTx = RegMethod.Execute(FixIm, MovIm)

which threw the error:

RuntimeError: 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(00000169FB1620A0): ComputeJacobianWithRespectToPosition not yet implemented for BSplineTransform

I found a thread about this error but it applied to ITK and didn’t seem directly applicable do to various calls that don’t seem to exist in SimpleITK .

I must be failing in understanding how to properly set up an initial BSpline transform using fiducials. I wonder if someone could please help explain what it is that I’m doing wrong or point me in the direction of an example.

For example, is the general idea to apply a “best fit” approach, e.g.

sitk.LandmarkBasedTransformInitializer(sitk.VersorRigid3DTransform(), FixPts, MovPts)

and composite transforms, or must I use an elastic approach, e.g.

LandmarkTx = sitk.LandmarkBasedTransformInitializerFilter()
...
LandmarkTx.SetBSplineNumberOfControlPoints(NumControlPts)

when initializing a BSpline transform with fiducials?

Thanks.

8

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?

9

Following @zivy’s help it turned out to be an issue with my choice in optimizer.

1 Like
10

I am trying to figure out how to use the bspline transformation for both image to image alignment and for alignment using landmarks. Both in 3D. Is there good documentation or examples from which I can learn how to do that? In particular, when using sitk.BSplineTransformInitializer(image1, transformDomainMeshSize, order), is the image necessary? I managed to work without an image for rigid and affine transformation.

Yoav

11

Hello @Yoav_Freund,

Welcome to SimpleITK!

An example illustrating the use of the BSplineTransform for 3D registration is available in the toolkit’s Jupyter notebook repository. The 6* notebook series deals with registration.

For the BSplineTransformInitializer we need the image because the BSplineTransform has a bounded spatial domain which in this case is defined by the whole image. The rigid/affine transformations have a global domain so we don’t need the image. For details on the various transformation types see this notebook.

12

Thanks Ziv!

I understand that to use splines you need to define a finite domain, but that seems different than specifying a whole 3D image. Can you just specify a bounding box, or do you need an actual 3D image?

I also understand that images are necessary when the optimization is based on similarity between images. However, I want to find a transformation that fits best a set of landmark pairs, so an image is not needed.

Do you think that using a 3d image with all zeros would work?

Finally, is there documentation that I can read, rather than example notebooks? The problem with the example notebooks is that they are very sparsely commented, so I am back to growing the official C++ API which I find very hard to use.

Yoav

13

Hi Yoav,

Yes, defining an arbitrary finite domain is different from using the whole image. Because the toolkit’s focus is on images, we implicitly assume that there is an image and that we are interested in its whole domain (the image’s physical extent).

To define an arbitrary domain, just create an appropriate image (defaults to all zeros) and feed that to the initializer:

x_size = 256
y_size = 128
z_size = 64
dummy = sitk.Image([x_size, y_size_, z_size], sitk.sitkUInt8)
dummy.SetOrigin([2,4,8])
dummy.SetSpacing((0.5,0.5,2))

Documentation that introduces the toolkit fundamental concepts, common conventions, installation instructions, and short examples can be found on read-the-docs. Hopefully these are easier to digest than the C++ docs or the notebooks.