Free form deformation or B-spline

1

I am trying to perform multi-resolution free form deformation/B-spline between two 2D images by following this: http://insightsoftwareconsortium.github.io/SimpleITK-Notebooks/Python_html/65_Registration_FFD.html

I don’t have any fixed or moving points as the images aren’t annotated.

outTx = bspline_intra_modal_registration(fixed_image=fixed_image,
                                      moving_image=out,
                                      # fixed_image_mask = (masks[fixed_image_index] == lung_label),
                                      # fixed_points = points[fixed_image_index],
                                      # moving_points = points[moving_image_index]
                                     )
resampler = sitk.ResampleImageFilter()
resampler.SetReferenceImage(fixed_image)
resampler.SetInterpolator(sitk.sitkNearestNeighbor)
resampler.SetDefaultPixelValue(0)
resampler.SetTransform(outTx)
out = resampler.Execute(out)

But I am not getting the 3 stage iteration like the example:
image
Also the fixed and moving image sizes are: fixed moving size: (66, 57) (57, 86)
On top the images don’t change at all.
Moving image before registration:
image
Moving image after registration:
image
The fixed image is:
image

Where am I making mistakes? I also changed the grid_physical_spacing but nothing significant happens.
MWE:

def command_iteration2(method):
    global metric_values
    print("metric value: {}".format(method.GetMetricValue()))
    metric_values.append(method.GetMetricValue())
metric_values = []
registration_method = sitk.ImageRegistrationMethod()
grid_physical_spacing = [.2, .2, .2]#, 50.0] # A control point every 50mm
image_physical_size = [size*spacing for size, spacing in zip(fixed_image.GetSize(), fixed_image.GetSpacing())]
mesh_size = [int(image_size/grid_spacing + 0.5)\
                 for image_size, grid_spacing in zip(image_physical_size, grid_physical_spacing)]
initial_transform = sitk.BSplineTransformInitializer(image1=fixed_image,
                                                     transformDomainMeshSize=mesh_size, order=3)
registration_method.SetMovingInitialTransform(compositeTx)
registration_method.SetInitialTransformAsBSpline(initial_transform,
                                                 inPlace=True,
                                                 scaleFactors=[1, 2, 4])

# Multi-resolution framework.
registration_method.SetShrinkFactorsPerLevel(shrinkFactors=[4, 2, 1])
registration_method.SetSmoothingSigmasPerLevel(smoothingSigmas=[4, 2, 0.01])
# registration_method.SmoothingSigmasAreSpecifiedInPhysicalUnitsOn()
registration_method.SetInterpolator(sitk.sitkNearestNeighbor)
registration_method.SetMetricAsANTSNeighborhoodCorrelation(16)
registration_method.MetricUseFixedImageGradientFilterOff()  # not sure what it does but saw it in example
registration_method.SetOptimizerScalesFromPhysicalShift()
# registration_method.SetMetricSamplingStrategy(registration_method.NONE)
registration_method.SetOptimizerAsLBFGS2(solutionAccuracy=1e-2,
                                         numberOfIterations=250,
                                         deltaConvergenceTolerance=1e-6)

registration_method.AddCommand(sitk.sitkIterationEvent, lambda: command_iteration2(registration_method))
outTx = registration_method.Execute(fixed_image, moving_image)
resampler = sitk.ResampleImageFilter()
resampler.SetReferenceImage(fixed_image)
resampler.SetInterpolator(sitk.sitkNearestNeighbor)
resampler.SetDefaultPixelValue(0)
resampler.SetTransform(outTx)
out = resampler.Execute(out)
plt.plot(metric_values)
plt.ylabel("Registration metric")
plt.xlabel("Iterations")
plt.show()
2

Hello @banikr,

In terms of pixels, these images are tiny to begin with (66x57, 57x86). In this scenario a multi-scale approach does not make much sense as these sizes are what is expected at the tip of the pyramid, not its base. I would recommend a single resolution registration.

Also, the mesh_size for the BSpline in this code is not appropriate. With BSpline registration you control the deformation via a grid of control points that are overlaid onto the image domain. The local deformation is controlled by these and the density of the grid points is much lower than the pixel grid. In the provided code it appears that for 66 pixels there are 300 control points, this does not make sense. You would probably want a control point every 10 pixels, so a grid of let’s say 7x6 for the 66x57 image.

1 Like
3

Hi @zivy,
though I managed to work around the multi-resolution framework on tiny images which may not necessary, but I can’t explain the ‘strange’ morph in results.
out is the moving image which comes after composite transform of rigid transformation: compositeTx

metric_values = []
grid_physical_spacing = [8, 6]  # , .2] # A control point every x, y mm
image_physical_size = [size * spacing for size, spacing in zip(fixed_image.GetSize(), fixed_image.GetSpacing())]
mesh_size = [int(image_size / grid_spacing + 0.5) \
             for image_size, grid_spacing in zip(image_physical_size, grid_physical_spacing)]
mesh_size = [int(sz / 4 + 0.5) for sz in mesh_size]
registration_method = sitk.ImageRegistrationMethod()
initial_transform = sitk.BSplineTransformInitializer(image1=fixed_image,
                                                     transformDomainMeshSize=mesh_size,
                                                     order=3)
registration_method.SetMovingInitialTransform(compositeTx)
registration_method.SetInitialTransformAsBSpline(initial_transform,
                                                 inPlace=True,
                                                 scaleFactors=[1, 2, 4])
# registration_method.SetMetricAsMeanSquares() # results are bad...very morphed. 
registration_method.SetMetricAsANTSNeighborhoodCorrelation(20)
registration_method.SetInitialTransform(initial_transform)
#         # RANDOM
registration_method.SetMetricSamplingStrategy(registration_method.NONE)
# Multi-resolution framework.
registration_method.SetShrinkFactorsPerLevel(shrinkFactors=[4, 2, 1])
registration_method.SetSmoothingSigmasPerLevel(smoothingSigmas=[2, 1, 0])
registration_method.SmoothingSigmasAreSpecifiedInPhysicalUnitsOn()
registration_method.SetInterpolator(sitk.sitkNearestNeighbor)
registration_method.SetOptimizerScalesFromPhysicalShift(
    # centralRegionRadius=4,
    # smallParameterVariation=0.001
)

registration_method.SetOptimizerAsGradientDescent(
    learningRate=.05,
    numberOfIterations=1000,
    convergenceMinimumValue=1e-5,
    convergenceWindowSize=150,
    estimateLearningRate=registration_method.EachIteration,
)

registration_method.AddCommand(sitk.sitkIterationEvent, lambda: command_iteration2(registration_method))

outTx = registration_method.Execute(fixed_image, moving_image)  # moving_image
print('Optimizer\'s stopping condition, {}'.format(registration_method.GetOptimizerStopConditionDescription()))
resampler = sitk.ResampleImageFilter()
resampler.SetReferenceImage(fixed_image)
resampler.SetInterpolator(sitk.sitkNearestNeighbor)
resampler.SetDefaultPixelValue(0)
resampler.SetTransform(outTx)
out = resampler.Execute(out)

plt.plot(metric_values)
plt.ylabel("Registration metric")
plt.xlabel("Iterations")
plt.show()

The out/moving image is a three label mask, I tried to normalize it between [0, 1] and register, which doesn’t change much.

reference/fixed image(normalized between 0-1):
image

out/moving image before B-Spline:
image
output:
image

Some of the conditions or prints to check:
Optimizer's stopping condition, GradientDescentOptimizerv4Template: Convergence checker passed at iteration 149.
metric curve for B-Spline:
image

Removing multi-resolution framework also seems to result in bad output:

image

Any suggestion where to improve from here?