ITKElastix: SetNumberOfThreads() segfaults with bspline parameter map

Engineering
1

When increasing the number of threads (> 1) of the itk.ElastixRegistrationMethod, the registration method segfaults when using the bspline default parameter map.
This already fails with the given example: bspline example

Code to change:

// Call registration function
elastix_object = itk.ElastixRegistrationMethod.New(fixed_image_bspline, moving_image_bspline)
elastix_object.SetParameterObject(parameter_object)
elastix_object.SetLogToConsole(True)
elastix_object.SetNumberOfThreads(8) # DOESNT WORK WITH BSPLINE
elastix_object.SetOutputDirectory("output")
elastix_object.UpdateLargestPossibleRegion()

result_image_bspline = elastix_object.GetOutput()
result_transform_parameters = elastix_object.GetTransformParameterObject()

//result_image_bspline, result_transform_parameters = itk.elastix_registration_method(
//    fixed_image_bspline, moving_image_bspline,
//    parameter_object=parameter_object,
//    log_to_console=True)

Even when lowering the MaximumNumberOfSamplingAttempts (as it suggests), it still segfaults.

Output:

WARNING: The parameter “FixedInternalImagePixelType”, requested at entry number 0, does not exist at all.
The default value “float” is used instead.
WARNING: The parameter “MovingInternalImagePixelType”, requested at entry number 0, does not exist at all.
The default value “float” is used instead.
Installing all components.
InstallingComponents was successful.

ELASTIX version: 5.1.0
Command line options from ElastixBase:
-out output/
-threads 8
WARNING: The parameter “UseDirectionCosines”, requested at entry number 0, does not exist at all.
The default value “true” is used instead.

WARNING: The option “UseDirectionCosines” was not found in your parameter file.
From elastix 4.8 it defaults to true!
This may change the behavior of your registrations considerably.

Command line options from TransformBase:
-t0 unspecified, so no initial transform used
WARNING: The parameter “BSplineTransformSplineOrder”, requested at entry number 0, does not exist at all.
The default value “3” is used instead.
WARNING: The parameter “UseCyclicTransform”, requested at entry number 0, does not exist at all.
The default value “false” is used instead.

Reading images…
Reading images took 0 ms.

WARNING: the fixed pyramid schedule is not fully specified!
A default pyramid schedule is used.
WARNING: the moving pyramid schedule is not fully specified!
A default pyramid schedule is used.
Initialization of all components (before registration) took: 0 ms.
Preparation of the image pyramids took: 4 ms.

Resolution: 0
WARNING: The parameter “ShowExactMetricValue”, requested at entry number 0, does not exist at all.
The default value “false” is used instead.
WARNING: The parameter “UseMultiThreadingForMetrics”, requested at entry number 0, does not exist at all.
The default value “true” is used instead.
WARNING: The parameter “ShowExactMetricValue”, requested at entry number 0, does not exist at all.
The default value “false” is used instead.
WARNING: The parameter “UseMultiThreadingForMetrics”, requested at entry number 0, does not exist at all.
The default value “true” is used instead.
Setting the fixed masks took: 0 ms.
Setting the moving masks took: 0 ms.
WARNING: The parameter “UseRelativeWeights”, requested at entry number 0, does not exist at all.
The default value “false” is used instead.
WARNING: The parameter “FixedImageBSplineInterpolationOrder”, requested at entry number 0, does not exist at all.
The default value “1” is used instead.
WARNING: The parameter “UseRandomSampleRegion”, requested at entry number 0, does not exist at all.
The default value “false” is used instead.
WARNING: The parameter “NumberOfHistogramBins”, requested at entry number 0, does not exist at all.
The default value “32” is used instead.
WARNING: The parameter “NumberOfFixedHistogramBins”, requested at entry number 0, does not exist at all.
The default value “32” is used instead.
WARNING: The parameter “NumberOfMovingHistogramBins”, requested at entry number 0, does not exist at all.
The default value “32” is used instead.
WARNING: The parameter “FixedLimitRangeRatio”, requested at entry number 0, does not exist at all.
The default value “0.01” is used instead.
WARNING: The parameter “MovingLimitRangeRatio”, requested at entry number 0, does not exist at all.
The default value “0.01” is used instead.
WARNING: The parameter “FixedKernelBSplineOrder”, requested at entry number 0, does not exist at all.
The default value “0” is used instead.
WARNING: The parameter “MovingKernelBSplineOrder”, requested at entry number 0, does not exist at all.
The default value “3” is used instead.
WARNING: The parameter “UseFastAndLowMemoryVersion”, requested at entry number 0, does not exist at all.
The default value “true” is used instead.
WARNING: The parameter “UseJacobianPreconditioning”, requested at entry number 0, does not exist at all.
The default value “false” is used instead.
WARNING: The parameter “FiniteDifferenceDerivative”, requested at entry number 0, does not exist at all.
The default value “false” is used instead.
WARNING: The parameter “NumberOfSamplesForSelfHessian”, requested at entry number 0, does not exist at all.
The default value “100000” is used instead.
WARNING: The parameter “SP_A”, requested at entry number 0, does not exist at all.
The default value “20” is used instead.

WARNING: You have set MaximumNumberOfSamplingAttempts to 8.
This functionality is known to cause problems (stack overflow) for large values.
If elastix stops or segfaults for no obvious reason, reduce this value.
You may select the RandomSparseMask image sampler to fix mask-related problems.

WARNING: The parameter “SigmoidInitialTime”, requested at entry number 0, does not exist at all.
The default value “0” is used instead.
WARNING: The parameter “MaxBandCovSize”, requested at entry number 0, does not exist at all.
The default value “192” is used instead.
WARNING: The parameter “NumberOfBandStructureSamples”, requested at entry number 0, does not exist at all.
The default value “10” is used instead.
WARNING: The parameter “UseAdaptiveStepSizes”, requested at entry number 0, does not exist at all.
The default value “true” is used instead.
WARNING: The parameter “UseConstantStep”, requested at entry number 0, does not exist at all.
The default value “false” is used instead.
WARNING: The parameter “MaximumStepLengthRatio”, requested at entry number 0, does not exist at all.
The default value “1” is used instead.
WARNING: The parameter “MaximumStepLength”, requested at entry number 0, does not exist at all.
The default value “1” is used instead.
WARNING: The parameter “NumberOfGradientMeasurements”, requested at entry number 0, does not exist at all.
The default value “0” is used instead.
WARNING: The parameter “NumberOfJacobianMeasurements”, requested at entry number 0, does not exist at all.
The default value “1000” is used instead.
WARNING: The parameter “SigmoidScaleFactor”, requested at entry number 0, does not exist at all.
The default value “0.1” is used instead.
Elastix initialization of all components (for this resolution) took: 0 ms.
Initialization of AdvancedMattesMutualInformation metric took: 1 ms.
Initialization of TransformBendingEnergy metric took: 0 ms.
Starting automatic parameter estimation for AdaptiveStochasticGradientDescent …
WARNING: The parameter “ASGDParameterEstimationMethod”, requested at entry number 0, does not exist at all.
The default value “Original” is used instead.
Computing JacobianTerms …
Computing the Jacobian terms took 0.005527s
NumberOfGradientMeasurements to estimate sigma_i: 3
Sampling gradients …
Segmentation fault (core dumped)

2

Full code to reproduce:

import itk
import numpy as np
import matplotlib.pyplot as plt

def image_generator(x1, x2, y1, y2, bspline=False):
    image = np.zeros([100, 100], np.float32)
    for x in range(x1, x2):
        for y in range(y1, y2):
            if bspline:
                y += x
                if x > 99 or y > 99:
                    pass
                else:
                    image[y, x] = 1
            else:
                image[y, x] = 1
    image = itk.image_view_from_array(image)
    return image

# Create test images
fixed_image_bspline = image_generator(25,65,25,65)
moving_image_bspline = image_generator(5,55,5,40, bspline=True)

# Import Default Parameter Map
parameter_object = itk.ParameterObject.New()
default_affine_parameter_map = parameter_object.GetDefaultParameterMap('affine',4)
default_affine_parameter_map['FinalBSplineInterpolationOrder'] = ['1']
#parameter_object.AddParameterMap(default_affine_parameter_map)
default_bspline_parameter_map = parameter_object.GetDefaultParameterMap('bspline',4)
default_bspline_parameter_map['FinalBSplineInterpolationOrder'] = ['1']
parameter_object.AddParameterMap(default_bspline_parameter_map)
#parameter_object.SetParameter("MaximumNumberOfSamplingAttempts", "2") 

# Call registration function
elastix_object = itk.ElastixRegistrationMethod.New(fixed_image_bspline, moving_image_bspline)
elastix_object.SetParameterObject(parameter_object)
elastix_object.SetLogToConsole(True)
elastix_object.SetNumberOfThreads(8) # DOESNT WORK WITH BSPLINE
elastix_object.SetOutputDirectory("output")
elastix_object.UpdateLargestPossibleRegion() # REQUIRED

result_image_bspline = elastix_object.GetOutput()
result_transform_parameters = elastix_object.GetTransformParameterObject()

#result_image_bspline, result_transform_parameters = itk.elastix_registration_method(
#    fixed_image_bspline, moving_image_bspline,
#    parameter_object=parameter_object,
#    log_to_console=True)

# Load Transformix Object
transformix_object = itk.TransformixFilter.New(moving_image_bspline)
transformix_object.SetTransformParameterObject(result_transform_parameters)

# Update object (required)
transformix_object.UpdateLargestPossibleRegion()

# Results of Transformation
result_image_transformix = transformix_object.GetOutput()

# Plot images
fig, axs = plt.subplots(1,4, sharey=True, figsize=[30,30])
plt.figsize=[100,100]
axs[0].imshow(result_image_bspline)
axs[0].set_title('Result', fontsize=30)
axs[1].imshow(fixed_image_bspline)
axs[1].set_title('Fixed', fontsize=30)
axs[2].imshow(moving_image_bspline)
axs[2].set_title('Moving', fontsize=30)
axs[3].imshow(result_image_transformix)
axs[3].set_title('Transformix', fontsize=30)
plt.show()
3

I did some further testing and yes the multithreading works with the elastix binary (version 5.1.0) through the commandline for the default bspline parameter map.
It doesnt work with the itk-elastix pip package (version 0.19.2).