# Using LandmarkBasedTransformInitializer for BSpline registrations

**URL:** <https://discourse.itk.org/t/using-landmarkbasedtransforminitializer-for-bspline-registrations/4096>\
**Category:** Algorithms\
**Tags:** simpleitk, registration, python\
**Created:** [May 4, 2021, 7:57pm UTC](https://discourse.itk.org/t/using-landmarkbasedtransforminitializer-for-bspline-registrations/4096 "2021-05-04T19:57:17Z")\
**Posts on this page:** 1\
**Showing post:** 7

<div class="post-metadata">

**Author:** ![ctorti](https://discourse.itk.org/user_avatar/discourse.itk.org/ctorti/32/1782_2.png) [@ctorti](https://discourse.itk.org/u/ctorti)\
**Post date:** [May 24, 2021, 12:58pm UTC](https://discourse.itk.org/t/using-landmarkbasedtransforminitializer-for-bspline-registrations/4096/7 "2021-05-24T12:58:47Z")

</div>

@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:

```auto
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)`:

```auto
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:

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

```

and

```auto
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):

```auto
MeshSize = [NumControlPts] * FixIm.GetDimension()

InitialTx = sitk.BSplineTransformInitializer(FixIm, MeshSize)

```

and combined both into a composite transform:

```auto
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:

```auto
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](https://discourse.itk.org/t/bspline-multiresolution-error-while-setting-transform/3562) 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.

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

```

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

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

```

when initializing a BSpline transform with fiducials?

Thanks.

---

_[View the full topic](https://discourse.itk.org/t/using-landmarkbasedtransforminitializer-for-bspline-registrations/4096)._