How can I perform 3d image registration using SimpleElastix from Python?

Question

I'm trying to register two 3d volumes. An attempt at this can found here. The code first generates two different volumes, both containing exactly one sphere of radius 4. I'm then trying to align them using the default translation parameter map. However, as can be seen in the last line (and from the plots if running locally), the result volume is not aligned with the fixed volume at all. When attempting the same procedure, in 2d this time, the resulting image does appear to be correctly aligned with the fixed image, as can be seen here. Am I using the SimpleElastix API incorrectly? I looked through the Github repo of SimpleElastix, but I could not find any examples of 3d image registration (at least not using volumes generated in Python and then converting them to ITK images).

Code from 3d example:

vol1 = np.zeros((50, 50, 50))
for x in range(vol.shape[0]):
    for y in range(vol.shape[1]):
        for z in range(vol.shape[2]):
            vol1[x, y, z] = np.linalg.norm(np.subtract([x, y, z], [5, 3, 2])) < 4

vol2 = np.zeros((50, 50, 50))
for x in range(vol.shape[0]):
    for y in range(vol.shape[1]):
        for z in range(vol.shape[2]):
            vol1[x, y, z] = np.linalg.norm(np.subtract([x, y, z], [20, 30, 10])) < 4

img_a = sitk.GetImageFromArray(vol1)
img_b = sitk.GetImageFromArray(vol2)

parameterMap = sitk.GetDefaultParameterMap('translation')

itk_filter = sitk.ElastixImageFilter()
itk_filter.LogToConsoleOn()
itk_filter.SetFixedImage(img_a)
itk_filter.SetMovingImage(img_b)
itk_filter.SetParameterMap(parameterMap)
itk_filter.Execute()

result_vol = sitk.GetArrayFromImage(itk_filter.GetResultImage())

np.max(np.abs(vol1 - result_vol))

Code from 2d example:

vol1 = np.zeros((50, 50))
for x in range(vol1.shape[0]):
    for y in range(vol1.shape[1]):
        vol1[x, y] = np.linalg.norm(np.subtract([x, y], [20, 20])) < 4

vol2 = np.zeros((50, 50))
for x in range(vol2.shape[0]):
    for y in range(vol2.shape[1]):
        vol2[x, y] = np.linalg.norm(np.subtract([x, y], [4, 5])) < 4

img_a = sitk.GetImageFromArray(vol1)
img_b = sitk.GetImageFromArray(vol2)

parameterMap = sitk.GetDefaultParameterMap('translation')

itk_filter = sitk.ElastixImageFilter()
itk_filter.LogToConsoleOn()
itk_filter.SetFixedImage(img_a)
itk_filter.SetMovingImage(img_b)
itk_filter.SetParameterMap(parameterMap)
itk_filter.Execute()

result_vol = sitk.GetArrayFromImage(itk_filter.GetResultImage())

np.max(np.abs(vol1 - result_vol))

Answer 1

The registration probably fails because the spheres do not overlap. You could try to smooth the images a bit more using e.g. sitk.DiscreteGaussian(). However, in this artificial example it is highly likely that simply moving the spheres closer to each other would result in a good registration.

Note it is difficult to register binary images as there is very little gradient information in these kinds of images: Gradients are non-zero at the border only, and zero everywhere else because of the regions are flat intensity-wise.