@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.998927TransformDomainMeshSize: [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.04629690.0277059 0.999615 -0.0012846
0.0463148 0 0.998927VersorRigid3DTransform (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.