MultiMetricMultiResolutionRegistration fails when SetNumberOfThreads is used

[Svdvoort]

When using MultiMetricMultiResolutionRegistration with the object-oriented interface and using SetNumberOfThreads, ITKElastix will throw a segmentation fault when setting more than one thread.

Code to test, this will result in a segmentation fault (based on the example here ).
Using python 3.7.7 and ITK-Elastix 0.14.1

import itk

# Import Images
fixed_image = itk.imread('data/CT_2D_head_fixed.mha', itk.F)
moving_image = itk.imread('data/CT_2D_head_moving.mha', itk.F)

parameter_object = itk.ParameterObject.New()
parameter_object.AddParameterFile('data/parameters_Bspline_Multimetric.txt')

elastix_object = itk.ElastixRegistrationMethod.New(fixed_image,moving_image)
elastix_object.SetParameterObject(parameter_object)

elastix_object.SetLogToConsole(False)

# You can set this to 1, but setting any number > 1 will result in a segmentation fault
elastix_object.SetNumberOfThreads(2)

elastix_object.UpdateLargestPossibleRegion()
result_image = elastix_object.GetOutput()

[thewtex]

Hi @Svdvoort ,

Thanks for the report and nice reproducible test.

I tried with Python 3.7, ITKElastix main as of today, and I could not reproduce the crash.

Could you please:

• itk-elastix 0.14.2
• Let use know more about your operating system and hard architecture

Thanks!

[Svdvoort]

Unfortunately, I'm not able to install itk-elastix 0.14.2.
I'm using Ubuntu 18.04, and itk-elastix 0.14.2 is not available for Ubuntu 18.04.
This is perhaps another issue, but considering that this is quite a breaking change, perhaps versioning it as 0.15.0 would be better.

I'll see if I can find a system that I can test version 0.14.2, but since that would be a different system in any case that is perhaps not ideal.

[thewtex]

@Svdvoort thanks for testing and the feedback.

It sounds like we have an issue similar to SuperElastix/elastix#714

CC @N-Dekker @stefanklein

[N-Dekker]

Thanks for reporting this issue @Svdvoort, but sorry, I haven't been able to reproduce this issue so far. I just tried the equivalent in C++, using just the elastix C++ library (without ITKElastix):

#include <itkElastixRegistrationMethod.h>
#include <gtest/gtest.h>

GTEST_TEST(itkElastixRegistrationMethod, Example07MultimetricMultiImageRegistration)
{
  using ImageType = itk::Image<float>;

  // Import Images
  const auto fixed_image = itk::ReadImage<ImageType>("data/CT_2D_head_fixed.mha");
  const auto moving_image = itk::ReadImage<ImageType>("data/CT_2D_head_moving.mha");

  const auto parameter_object = elx::ParameterObject::New(); 
  parameter_object->AddParameterFile("data/parameters_Bspline_Multimetric.txt");
  
  const auto elastix_object = itk::ElastixRegistrationMethod<ImageType, ImageType>::New();
  elastix_object->SetFixedImage(fixed_image);
  elastix_object->SetMovingImage(moving_image);
  elastix_object->SetParameterObject(parameter_object);

  elastix_object->SetLogToConsole(false);

  // You can set this to 1, but setting any number > 1 will result in a segmentation fault
  elastix_object->SetNumberOfThreads(2);

  elastix_object->UpdateLargestPossibleRegion();
  const auto result_image = elastix_object->GetOutput();
}

Would you expect the segmentation fault to be reproducible by the C++ code above here? Honestly I only tried Windows/VS2019 Debug so far, but it seems to run fine.

[Svdvoort]

My C++ knowledge is a bit more limited, but from what I can tell this should indeed cause an issue.

Could it be because of older libraries? The new Elastix and ITK-Elastix don't run on Ubuntu 18.04 anymore due to outdated C libraries, so perhaps that's the underlying issue? It might therefore be resolved on Ubuntu 20.04/ with the latest ITK-elastix (0.14.2) but I'm not able to test that currently.

[Svdvoort]

To give an update on this:

I tried with itk-elastix 0.14.3 (this version can be installed on Ubuntu 18.04). Now the above code gives me errors both when setting 1 thread, and when setting multiple threads (in all cases a Segmentation Fault, without any additional output).

Therefore I also tried running with itk-elastix 0.14.1, but this gives an AttributeError: PyCapsule_Import "_ITKCommonPython._C_API" is not valid error. Similar to #154 and #159. However, the solution posted in #159 doesn't work anymore, because it seems that there is no 5.3rc4 version for itk-numerics anymore on pypi, only the post-releases. So I'm not currently able to reproduce the results as initially posted, but the latest version of itk-elastix gives me an error in all cases.

This is on Ubuntu 18.04, with python 3.7.7 in a virtual environment, with itk-elastix installed through pip 22.3.1

edit: This might be an issue with itk-elastix 0.14.3. It seems that no matter what code I run using itk-elastix (even just loading a parameter map) it results in a segmentation fault. Using an older environment in which itk-elastix 0.14.1 is still available works fine. I've made a new issue about this #177

[Svdvoort]

I've updated to itk-elastix 0.14.4, and this version is working again. Although the problem is not completely solved yet, it does output more than previously.

Before after running the above script it instantly gave an error: Segmentation fault (core dumped)

Now it does display a lot more output, but still ends up in the same error. The last few lines of output:

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.005726s
  NumberOfGradientMeasurements to estimate sigma_i: 2
  Sampling gradients ...
Segmentation fault (core dumped)

So it seems like the error occurs in the sampling of the gradients.

[stephahn0916]

Have the same behavior with itk-elastix 0.14.4 on macos.

[Svdvoort]

Issue still persists with itk-elastix 0.15.0, with the same output as version 0.14.4.

[Svdvoort]

To update on this status and hopefully focus in on the problem, here some extra information and tests. The problem still seems to persist, even when using the new MultiThreaderBase. The threading does seem to work properly with "native" itk functions

Test environment

The following examples and results have all been tested on Ubuntu 22.04.4, using Python 3.11.9. This has been tested with two versions of itk-elastix:

1. Through pypi with pip install itk-elastix
2. Downloading the wheel LinuxWheel311_2_28-x64 from the latest action at the time of writing

In both cases the results were the same, so all results will only be shown once. Below the packages that are installed in the Python environment:

itk==5.4.0
itk-core==5.4.0
itk-filtering==5.4.0
itk-io==5.4.0
itk-numerics==5.4.0
itk-registration==5.4.0
itk-segmentation==5.4.0
numpy==2.0.0

In the case of the pip install this also included itk-elastix==0.20.0, in the case of the wheel install it included: itk-elastix @ file:///LinuxWheel311_2_28-x64/itk_elastix-0.20.0-cp311-abi3-manylinux_2_28_x86_64.whl#sha256=4e910fcbbf86db13909471da1878b7e40385dfb94111a57e5b38270caafff559

Elastix tests

SetNumberOfThreads

The following script still gives a segmentation fault:

import itk

# Import Images
fixed_image = itk.imread('data/CT_2D_head_fixed.mha', itk.F)
moving_image = itk.imread('data/CT_2D_head_moving.mha', itk.F)

parameter_object = itk.ParameterObject.New()
parameter_object.AddParameterFile('data/parameters_Bspline_Multimetric.txt')

elastix_object = itk.ElastixRegistrationMethod.New(fixed_image,moving_image)
elastix_object.SetParameterObject(parameter_object)

elastix_object.SetLogToConsole(True)

# You can set this to 1, but setting any number > 1 will result in a segmentation fault
elastix_object.SetNumberOfThreads(2)

elastix_object.UpdateLargestPossibleRegion()
result_image = elastix_object.GetOutput()

This is the same behaviour as before and is therefore unchanged.

MultiThreaderBase

SetGlobalDefaultNumberOfThreads

Example script to work with the new MultiThreaderBase:

import itk

threader = itk.MultiThreaderBase.New()
threader.SetGlobalDefaultNumberOfThreads(2)

# Import Images
fixed_image = itk.imread("data/CT_2D_head_fixed.mha", itk.F)
moving_image = itk.imread("data/CT_2D_head_moving.mha", itk.F)

parameter_object = itk.ParameterObject.New()
parameter_object.AddParameterFile("data/parameters_Bspline_Multimetric.txt")

elastix_object = itk.ElastixRegistrationMethod.New(fixed_image, moving_image)
elastix_object.SetParameterObject(parameter_object)

elastix_object.SetLogToConsole(True)

filter_default_threads = (
    elastix_object.GetMultiThreader().GetGlobalDefaultNumberOfThreads()
)
print(filter_default_threads)

elastix_object.UpdateLargestPossibleRegion()
result_image = elastix_object.GetOutput()

This new option works and the print does show that the number of thread is set to 2. However, something doesn't seem to work properly. When using threader.SetGlobalDefaultNumberOfThreads(1) and monitoring resources with htop, it nicely shows only one process running. However setting anything > 1 shows that all cores/threads available are being used.

SetGlobalMaximumNumberOfThreads

Below test script using SetGlobalMaximumNumberOfThreads:

import itk

threader = itk.MultiThreaderBase.New()
threader.SetGlobalMaximumNumberOfThreads(2)

# Import Images
fixed_image = itk.imread("data/CT_2D_head_fixed.mha", itk.F)
moving_image = itk.imread("data/CT_2D_head_moving.mha", itk.F)

parameter_object = itk.ParameterObject.New()
parameter_object.AddParameterFile("data/parameters_Bspline_Multimetric.txt")

elastix_object = itk.ElastixRegistrationMethod.New(fixed_image, moving_image)
elastix_object.SetParameterObject(parameter_object)

elastix_object.SetLogToConsole(True)

elastix_object.UpdateLargestPossibleRegion()
result_image = elastix_object.GetOutput()

This gives the same result as SetGlobalDefaultNumberOfThreads: more than the specified 2 threads are used.

itk tests

The following example code uses an itk filter, and here setting the number of threads through MultiThreaderBase does work properly. Setting the number of threads to 2 shows in htop that only 2 threads are being used, also verified this for setting it to 1 and 4 threads. This is the behaviour I would expect, but somehow does not seem to generalize to the elastix functions. Using SetGlobalMaximumNumberOfThreads instead of SetGlobalDefaultNumberOfThreads in this example results in the same behaviour.

import itk

threader = itk.MultiThreaderBase.New()
threader.SetGlobalDefaultNumberOfThreads(2)


dimension = 3
float_image_type = itk.Image[itk.F, dimension]


corner = itk.Index[dimension]()
corner.Fill(0)

size = itk.Size[dimension]()
size.Fill(500)

region = itk.ImageRegion[dimension]()
region.SetIndex(corner)
region.SetSize(size)

image = float_image_type.New(Regions=region)
image.Allocate()
image.FillBuffer(0)

# Make a square
image[40:100, 40:100] = 100

# Define some types
unsigned_char_image_type = itk.Image[itk.UC, dimension]
filt = itk.MedianImageFilter[float_image_type, float_image_type].New()
filt.SetInput(image)
filt.SetRadius(3)
filt.Update()

Not initiating MultiThreaderBase

SetGlobalDefaultNumberOfThreads

There is some interesting behaviour however when not calling itk.MultiThreaderBase.New() first as in the following script:

import itk

threader = itk.MultiThreaderBase.SetGlobalDefaultNumberOfThreads(2)

dimension = 3
float_image_type = itk.Image[itk.F, dimension]


corner = itk.Index[dimension]()
corner.Fill(0)

size = itk.Size[dimension]()
size.Fill(500)

region = itk.ImageRegion[dimension]()
region.SetIndex(corner)
region.SetSize(size)

image = float_image_type.New(Regions=region)
image.Allocate()
image.FillBuffer(0)

# Make a square
image[40:100, 40:100] = 100

# Define some types
unsigned_char_image_type = itk.Image[itk.UC, dimension]
filt = itk.MedianImageFilter[float_image_type, float_image_type].New()
filt.SetInput(image)
filt.SetRadius(3)
filt.Update()

SetGlobalMaximumNumberOfThreads

The above script works as expected, however using SetGlobalMaximumNumberOfThreads results in the same behaviour as for the elastix-filter: all available cores/threads are being used. The following is the test script:

import itk

threader = itk.MultiThreaderBase.SetGlobalMaximumNumberOfThreads(2)

dimension = 3
float_image_type = itk.Image[itk.F, dimension]


corner = itk.Index[dimension]()
corner.Fill(0)

size = itk.Size[dimension]()
size.Fill(500)

region = itk.ImageRegion[dimension]()
region.SetIndex(corner)
region.SetSize(size)

image = float_image_type.New(Regions=region)
image.Allocate()
image.FillBuffer(0)

# Make a square
image[40:100, 40:100] = 100

# Define some types
unsigned_char_image_type = itk.Image[itk.UC, dimension]
filt = itk.MedianImageFilter[float_image_type, float_image_type].New()
filt.SetInput(image)
filt.SetRadius(3)
filt.Update()