Different BSpline fixed-parameters for equal-sized images?

Beginner Questions
,
1

I started with a very simple registration pipeline, basically like this:

import itk

fixed_image = itk.imread('fixed.mhd', itk.F)
moving_image = itk.imread('moving.mhd', itk.F)

fixed_mask = itk.imread("fixed_mask.mhd", itk.UC)
moving_mask = itk.imread("moving_mask.mhd", itk.UC)

RESOLUTIONS = [4, 2, 1]

parameter_object = itk.ParameterObject.New()

rigid_map = itk.ParameterObject.GetDefaultParameterMap("rigid", len(RESOLUTIONS))
rigid_map['Transform'] = ["SimilarityTransform"]
rigid_map["AutomaticTransformInitialization"] = ["true"]
rigid_map["AutomaticTransformInitializationMethod"] = ["GeometricalCenter"]
rigid_map["Metric"] =  ["AdvancedNormalizedCorrelation"]
parameter_object.AddParameterMap(rigid_map)

spline_map = itk.ParameterObject.GetDefaultParameterMap("bspline", len(RESOLUTIONS), 20.0)
spline_map["Metric"] =  ('AdvancedNormalizedCorrelation', 'TransformBendingEnergyPenalty')
parameter_object.AddParameterMap(spline_map)
parameter_object.SetParameter("ImagePyramidSchedule", [str(x) for x in RESOLUTIONS for _ in range(3)])

# Load Elastix Image Filter Object
registration_method= itk.ElastixRegistrationMethod.New(
    fixed_image=fixed_image,
    moving_image=moving_image,
    parameter_object=parameter_object,
    fixed_mask=fixed_mask,
    moving_mask=moving_mask,
)
registration_method.UpdateLargestPossibleRegion()
itk_composite_transform = registration_method.ConvertToItkTransform(registration_method.GetCombinationTransform())

so far this works great. However, I found out that even if my fixed and moving images all have the same sizes, I get different BSpline grid origins and spacings.
For example, I’m registering 19 different images onto one, then these are my fixed parameters:

[5.0, 5.0, 8.0, -25.35798058370918, -26.53888280642838, -14.484514672372548, 20.39789151360511, 19.57830637393795, 18.307260209617585, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0]
[5.0, 5.0, 8.0, -24.90272869368449, -29.97265966883625, -14.597313349700325, 19.870291464945765, 21.007113881112293, 18.147430027811183, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0]
[5.0, 5.0, 8.0, -27.279253675326395, -47.45167550654084, -21.293408546504715, 21.290357264138954, 30.386644116875424, 20.510496013529377, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0]
[5.0, 5.0, 8.0, -24.730771152294917, -26.875317776098214, -15.907550783130628, 20.117265133642313, 19.600241501980335, 18.667774405055834, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0]
[5.0, 5.0, 8.0, -31.86043813671492, -32.9392057081962, -20.36568622166342, 23.42378658667599, 22.44038606652115, 19.966109186614773, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0]
[5.0, 5.0, 8.0, -23.340586513524098, -26.99395602983922, -18.32371078159948, 19.65497064383643, 20.018114832493765, 19.39686744650334, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0]
[5.0, 5.0, 8.0, -23.27130007768004, -28.70194159153989, -17.14189213754565, 19.488016541023367, 20.38525182520469, 18.901488866675557, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0]
[5.0, 5.0, 8.0, -23.06263411579611, -25.52238763479419, -13.016831616922719, 19.07277330143413, 18.866759709815124, 18.342552858211178, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0]
[5.0, 5.0, 8.0, -23.222969896798983, -25.43296759085781, -14.808330793497865, 19.31795475279227, 18.962923357017615, 18.316815929901797, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0]
[5.0, 5.0, 8.0, -25.22324817980321, -27.03354767851697, -14.59571408076085, 20.516447936958066, 19.623321712998102, 18.46156198360738, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0]
[5.0, 5.0, 8.0, -31.51386642669825, -34.003003909507505, -21.657093792441675, 22.70194405940115, 23.22194074637386, 20.244822158617023, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0]
[5.0, 5.0, 8.0, -27.64451431272127, -26.72970564291126, -19.262585468839767, 21.05932734294045, 19.68177490483581, 19.409880913478865, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0]
[5.0, 5.0, 8.0, -24.848542480606426, -33.91185783248196, -17.67613442069956, 19.88833664938058, 22.592449607854608, 18.954712769421022, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0]
[5.0, 5.0, 8.0, -30.05597062852287, -30.87491702570174, -18.10566480841071, 21.710463607575576, 22.101068159946635, 19.235312235255364, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0]
[5.0, 5.0, 8.0, -34.87718554047299, -31.926563782949774, -22.87149339867384, 23.913435003312465, 22.69052850596975, 20.748620132845936, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0]
[5.0, 5.0, 8.0, -24.04544329389816, -27.85637307997677, -16.45823455204542, 19.395408267458137, 19.85101796649231, 18.854031755866696, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0]
[5.0, 5.0, 8.0, -24.055699774968087, -27.03061740016157, -19.224972452579408, 19.43841874219105, 19.84148900872092, 19.340374043121024, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0]
[5.0, 5.0, 8.0, -26.640964862302116, -29.91395913673105, -20.5672273913215, 20.701147414389467, 20.54574679419675, 19.9045059736219, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0]
[5.0, 5.0, 8.0, -23.772821471356334, -29.618726892156687, -17.278169788390386, 19.92443589543942, 20.83295970787472, 19.109609329066828, 1.0, 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 1.0]

As you can see, they all have the same number of knot spans (and direction matrix), but their grid origins and grid spacings vary.
Is this maybe because of the masks?
Or is the issue that I use a SimilarityTransform? (btw, is it correct that this does not have a default parameter map?) The documentation says something about affine transforms:

If an affine transformation is provided as initial transformation, the control grid will be scaled to cover the fixed image domain in the space defined by the initial transformation.

However, I’m also not providing an initial transform…

Is there any parameter to fix those, so for an equally-sized image, I get the same grid?

2

So, I tested now on some synthetic images. Interstingly, if I only add a bspline to the registration parameter map, I get always the same fixed parameters. It is only if I add a rigid transform in the beginning, that the parameters change. So is this the hint from the documentation about initial transforms? Is there any way to change this, such that all my splines have the same fixed parameters?

Here is a MWE:

import itk
import numpy as np

x = np.zeros((100, 100, 100))
x[20:80, 20:80, 20:80] = 1

X, Y, Z = np.mgrid[0:100, 0:100, 0:100]
# Add some "random" pattern...
x *= np.cos(X / 4) + np.sin(Y / 4) + np.cos(Z / 4) * np.sin(Z / 4)
# Add a border
x[18:82, 18:82, 18:82] += 2
x[20:80, 20:80, 20:80] -= 2
fixed = itk.GetImageFromArray(np.ascontiguousarray(x))


def get_random_bspline(img):
    transform = itk.BSplineTransform[itk.D, 3, 3].New()
    initializer = itk.BSplineTransformInitializer[BTransformType, type(img)].New()
    initializer.SetTransform(transform)
    initializer.SetImage(img)
    initializer.SetTransformDomainMeshSize([5, 5, 5])
    initializer.InitializeTransform()
    transform.SetIdentity()
    
    parameters = transform.GetParameters()
    for i in range(transform.GetNumberOfParameters()):
        parameters[i] = float(np.random.normal(0, 10))
    transform.SetParameters(parameters)

    return transform


moving = [
    itk.resample_image_filter(
        fixed,
        interpolator=itk.LinearInterpolateImageFunction.New(fixed),
        transform=get_random_bspline(fixed),
        reference_image = fixed,
        use_reference_image = True,
    )
    for _ in range(10)
]

for moving_image in moving:   
    parameter_object = itk.ParameterObject.New()
    # --> This seems to be causing a different grid setting...
    # parameter_object.AddParameterMap(itk.ParameterObject.GetDefaultParameterMap("rigid", 3))
    parameter_object.AddParameterMap(itk.ParameterObject.GetDefaultParameterMap("bspline", 3, 25.0))
    
    registration_method= itk.ElastixRegistrationMethod.New(
        fixed_image=fixed,
        moving_image=moving_image,
        parameter_object=parameter_object,
    )
    registration_method.UpdateLargestPossibleRegion()
    tfm = itk.CompositeTransform[itk.D, 3].cast(
        registration_method.ConvertToItkTransform(registration_method.GetCombinationTransform())
    )

    print(*map('{:.2f}'.format, tfm.GetNthTransform(0).GetFixedParameters()), sep=', ')

Without the rigid transform, I get:

7.00, 7.00, 7.00, -25.50, -25.50, -25.50, 25.00, 25.00, 25.00, 1.00, 0.00, 0.00, 0.00, 1.00, 0.00, 0.00, 0.00, 1.00
7.00, 7.00, 7.00, -25.50, -25.50, -25.50, 25.00, 25.00, 25.00, 1.00, 0.00, 0.00, 0.00, 1.00, 0.00, 0.00, 0.00, 1.00
7.00, 7.00, 7.00, -25.50, -25.50, -25.50, 25.00, 25.00, 25.00, 1.00, 0.00, 0.00, 0.00, 1.00, 0.00, 0.00, 0.00, 1.00
7.00, 7.00, 7.00, -25.50, -25.50, -25.50, 25.00, 25.00, 25.00, 1.00, 0.00, 0.00, 0.00, 1.00, 0.00, 0.00, 0.00, 1.00
7.00, 7.00, 7.00, -25.50, -25.50, -25.50, 25.00, 25.00, 25.00, 1.00, 0.00, 0.00, 0.00, 1.00, 0.00, 0.00, 0.00, 1.00
7.00, 7.00, 7.00, -25.50, -25.50, -25.50, 25.00, 25.00, 25.00, 1.00, 0.00, 0.00, 0.00, 1.00, 0.00, 0.00, 0.00, 1.00
7.00, 7.00, 7.00, -25.50, -25.50, -25.50, 25.00, 25.00, 25.00, 1.00, 0.00, 0.00, 0.00, 1.00, 0.00, 0.00, 0.00, 1.00
7.00, 7.00, 7.00, -25.50, -25.50, -25.50, 25.00, 25.00, 25.00, 1.00, 0.00, 0.00, 0.00, 1.00, 0.00, 0.00, 0.00, 1.00
7.00, 7.00, 7.00, -25.50, -25.50, -25.50, 25.00, 25.00, 25.00, 1.00, 0.00, 0.00, 0.00, 1.00, 0.00, 0.00, 0.00, 1.00
7.00, 7.00, 7.00, -25.50, -25.50, -25.50, 25.00, 25.00, 25.00, 1.00, 0.00, 0.00, 0.00, 1.00, 0.00, 0.00, 0.00, 1.00

So all have the same grid, and as you can see the grid is always 25mm - as specified.
However, when the rigid is added, I get:

7.00, 7.00, 7.00, -32.83, -31.69, -31.29, 27.12, 26.44, 26.98, 1.00, 0.00, 0.00, 0.00, 1.00, 0.00, 0.00, 0.00, 1.00
7.00, 7.00, 7.00, -24.08, -25.27, -23.88, 25.22, 25.08, 25.27, 1.00, 0.00, 0.00, 0.00, 1.00, 0.00, 0.00, 0.00, 1.00
7.00, 7.00, 7.00, -29.01, -27.00, -30.76, 26.67, 25.70, 27.31, 1.00, 0.00, 0.00, 0.00, 1.00, 0.00, 0.00, 0.00, 1.00
7.00, 7.00, 7.00, -28.35, -28.82, -33.14, 26.09, 26.25, 27.19, 1.00, 0.00, 0.00, 0.00, 1.00, 0.00, 0.00, 0.00, 1.00
7.00, 7.00, 7.00, -26.57, -30.02, -30.12, 25.21, 26.47, 26.27, 1.00, 0.00, 0.00, 0.00, 1.00, 0.00, 0.00, 0.00, 1.00
7.00, 7.00, 7.00, -29.79, -32.08, -29.60, 26.18, 26.40, 26.17, 1.00, 0.00, 0.00, 0.00, 1.00, 0.00, 0.00, 0.00, 1.00
7.00, 7.00, 7.00, -32.12, -30.06, -29.41, 27.01, 27.16, 26.48, 1.00, 0.00, 0.00, 0.00, 1.00, 0.00, 0.00, 0.00, 1.00
7.00, 7.00, 7.00, -26.17, -26.78, -25.33, 26.21, 25.99, 25.34, 1.00, 0.00, 0.00, 0.00, 1.00, 0.00, 0.00, 0.00, 1.00
7.00, 7.00, 7.00, -31.43, -29.89, -30.68, 26.68, 26.36, 26.75, 1.00, 0.00, 0.00, 0.00, 1.00, 0.00, 0.00, 0.00, 1.00
7.00, 7.00, 7.00, -31.21, -27.75, -29.61, 26.80, 26.33, 26.31, 1.00, 0.00, 0.00, 0.00, 1.00, 0.00, 0.00, 0.00, 1.00

Which I find a bit weird, as the moving and fixed domain are exactly the same - thus why does it need to adjust the bspline grid?

Edit: Okay, so the code that is called here seems to be elastix/Common/Transforms/itkGridScheduleComputer.hxx at 2aafd7b0a4e73259523c2d021b92edb37df2c890 · SuperElastix/elastix · GitHub
But, I do not see any setting to turn this off - except by not providing a initial transform. So the only way would be to run the rigid registration, transform the image and run the bspline registration on them?

3

Will this not also produce differing BSpline parameters, as your resampled images will have slightly different origins and maybe direction matrices?

A possibility is to use the same rigid transform for all your images, if that would work. If the rigid transform is not radical, maybe you can just skip it and apply BSpline only.

4

I may have misunderstood something in registration fundamentally wrong in the last years that I use it, but if the fixed domain and the moving domain are the same (i.e., all images have same origin, size, spacing and direction), then why should the resampled images have different origins or directions? The moving image should always be resampled to the fixed domain, no?
Of course, it could happen that the object moves outside the image domain due to the affine transform - which I guess is the reason for this adaptation of the bspline domain?

I fear that would not work. At least some translation and isotropic scaling is often needed to produce proper overlap…

5

With rigid and affine transforms, images do not need to be resampled. Take a look at TransformGeometryImageFilter. I guess that elastix does this, as this transformation is lossless. There is not resampling loss, and nothing needs to be cut off, as nothing goes beyond the edge.

If you resample images to the same image grid after rigid/affine transform, then do deformable registration as a separate step, your BSpline definition coefficients should all be the same.

6

aha :slight_smile: I always learn new tricks with ITK :wink:
thanks, this totally makes sense now!

yep, okay, this was the procedure I had implemented so far, and will adopt for elastix too.

1 Like
7

So I implemented it now as a two stage process with in-between resampling, which works as expected!
However, if there is another solution or parameter I may have overlooked, I would be happy to know about it :wink:

1 Like