ENH: Make ElastixMain database creation + loading components thread safe

[N-Dekker]

Ensures that only one thread may load the components into the ComponentDatabase, using using C++11 "magic statics".

ElastixMain::GetComponentDatabase() now returns a "const reference" instead of a non-const pointer, to prevent one thread from modifying the database while another is reading from the database.

Removed the two UnloadComponents() member functions (from both ElastixMain and ComponentLoader), as they appear not useful anymore.

Aims to fix issue #174
"elastix 4.9 static libary version not thread safe"
reported by @jiangliMED, August 16, 2019

With help from Konstantinos Ntatsis (@ntatsisk)

[N-Dekker]

@ntatsisk Hi Konstantinos! Can you please review this pull request? And could you possibly also try it out, to see if it does indeed fix a multi-threading issue that you encountered? You may check out the (temporary) branch of this pull request: https://github.com/SuperElastix/elastix/tree/ElastixMain-LoadComponents-thread-safe

[ntatsisk]

Hello @N-Dekker, the solution looks great! I really liked the combination of magic statics together with a lambda function for the ComponentDatabase initialization. I merged this branch with my local branch where the logger is "deactivated" and managed to run some tests. It works with no issues!

[N-Dekker]

@ntatsisk Thank you for your approval and your enthusiastic reply! I certainly enjoyed combining magic statics together with a lambda 😃 May I ask what your use case looks like? So you call RegisterImages multiple times simultaneously within the same process? How do you do that? Something like std::async? And what do you gain from doing so?

[ntatsisk]

@N-Dekker No, we are using OpenMP to parallelize the registrations. The idea is as follows:

#pragma omp parallel for default(none) num_threads(24)
	for (int i = 0; i < 600; i++)
	{
			ElastixFilterType::Pointer elastixFilter = ElastixFilterType::New();
			elastixFilter->SetFixedImage(images[i]);
			elastixFilter->SetMovingImage(images[i+1]);
			elastixFilter->SetParameterObject(m_elxParameterObjectPointer);
			elastixFilter->LogToConsoleOff();
			elastixFilter->LogToFileOff();

// #pragma omp critical // (= not necessary if elastix is thread-safe)
			elastixFilter->Update();
	}

As you see, we are trying to do many (simple) registrations but in a short amount of time (~<10sec) which is a requirement. So, it is better to split the registrations among the threads instead of using many threads for a single registration, as ITK threading does. The gain is some saved seconds eg. from 12sec (with the critical section) to 7sec (in the thread-safe version). In another part, it goes from 30sec to 10sec. Generally, these times might not seem like a big deal, but they are for our application ;)

[N-Dekker]

@mstaring @stefanklein Hope it's OK to you that I'll merge this pull request this afternoon!

[mstaring]

I do not get all the details of this PR, but i assume i) elastix still works (all tests pass for the command line version and ITKElastix also still works), and ii) it was verified that it actually fixes a problem (a new test first failed but now passes).

[N-Dekker]

@mstaring Thanks for your feedback

i assume i) elastix still works (all tests pass for the command line version and ITKElastix also still works)

All elastix test are passed successfully, at the CI. However, I did not test ITKElastix with this pull request. We currently do not test ITKElastix on each elastix pull request.

But I don't think ITKElastix would have any problem with this pull request. Unless it would try to call UnloadComponents(), which is removed by this pull request. But that seems unlikely to me, especially because UnloadComponents() wasn't useful anymore anyway. @ViktorvdValk do you know if UnloadComponents() was wrapped to Python, by ITKElastix?

ii) it was verified that it actually fixes a problem (a new test first failed but now passes).

It is verified by Konstantinos (@ntatsisk) that it actually fixes a problem. I can still see if an automatic test can be added, based on his use case.

[mstaring]

sounds good then!

[N-Dekker]

@ntatsisk Sorry I did not yet merge, because I'm still trying to write a unit test (GoogleTest) for this pull request. But if you need it now, I can just merge this and write the test later. What do you think? Anyway, I'm just trying (locally on my laptop):

GTEST_TEST(ElastixFilter, Parallel)
{
  using ImageType = itk::Image<unsigned char>;
  std::vector<ImageType::Pointer> images(601);

  for (auto & image : images)
  {
    image = ImageType::New();
    image->SetRegions({ 100, 100 });
    image->Allocate(true);
  }

#pragma omp parallel for
  for (int i{}; i < 600; ++i)
  {
    const auto elastixFilter = elx::ElastixFilter<ImageType, ImageType>::New();
    elastixFilter->SetFixedImage(images[i]);
    elastixFilter->SetMovingImage(images[i + 1]);
    elastixFilter->LogToConsoleOff();
    elastixFilter->LogToFileOff();
    elastixFilter->Update();
  }
}

But it says:

User Error 1001: omp_set_num_threads should only be called in serial regions

Do you have a suggestion?

[ntatsisk]

Hi @N-Dekker, I am not in a hurry at all with the merging so it can wait. Regarding the error, I assume you are compiling with "warnings as errors" because I have got the same message but as a warning instead. I used to get it while running in Debug mode but it has been a while since then. Maybe, it is coming from here:

elastix/Components/Optimizers/StandardGradientDescent/itkGradientDescentOptimizer2.cxx

Lines 212 to 224 in af382ee

	#else // Otherwise use OpenMP
	/** Get a reference to the current position. */
	const ParametersType & currentPosition = this->GetScaledCurrentPosition();
	/** Update the new position. */
	const int nthreads = static_cast<int>(this->m_Threader->GetNumberOfWorkUnits());
	omp_set_num_threads(nthreads);
	# pragma omp parallel for
	for (int j = 0; j < static_cast<int>(spaceDimension); j++)
	{
	newPosition[j] = currentPosition[j] - this->m_LearningRate * this->m_Gradient[j];
	}
	#endif

Since this omp_set_num_threads is now called inside the user's omp parallel for directive (like the one that you have in your test script), the warning is emitted. To be honest, I haven't fully understood the internal elastix threading and whether this and the other few OpenMP clauses are necessary/helpful. If they do need to stay then I would try re-writing the clause in the above link as follows:

#pragma omp parallel for num_threads(nthreads)

Does this work? Btw, even if it does, I think your test will still crash due to the logger being thread-unsafe, right?

[N-Dekker]

Regarding the error, I assume you are compiling with "warnings as errors" because I have got the same message but as a warning instead. I used to get it while running in Debug mode but it has been a while since then. Maybe, it is coming from here:

elastix/Components/Optimizers/StandardGradientDescent/itkGradientDescentOptimizer2.cxx

Lines 212 to 224 in af382ee

	#else // Otherwise use OpenMP
	/** Get a reference to the current position. */
	const ParametersType & currentPosition = this->GetScaledCurrentPosition();
	/** Update the new position. */
	const int nthreads = static_cast<int>(this->m_Threader->GetNumberOfWorkUnits());
	omp_set_num_threads(nthreads);
	# pragma omp parallel for
	for (int j = 0; j < static_cast<int>(spaceDimension); j++)
	{
	newPosition[j] = currentPosition[j] - this->m_LearningRate * this->m_Gradient[j];
	}
	#endif

Thanks for the suggestion, Konstantinos, but as far as I can see, those lines of code are excluded from compilation, as they are part of the "else block" of an #if 1!

Anyway, I just submitted a new "ElastixFilter.UpdateInParallel" test: #393 My aim is to have a little test that just shows that this pull request (#389) really fixes a certain problem. But the test is still unstable! I see it crashing some times. Could you also please have a look?

[ntatsisk]

Ops, my bad! So, it is clearly not the point that is called (maybe this instead?). Anyway, I also noticed that, actually, all the #pragma omp parallel directives are manually deactivated from the entire elastix codebase at this moment.

Regarding the crashing, as I mentioned in my last comment, I think it is because of the logger ;)

[N-Dekker]

Regarding the crashing, as I mentioned in my last comment, I think it is because of the logger ;)

I believe so too. The logger is still troublesome. So then, is there still any possible test case for this specific pull request, which only improves the thread-safety of the component database?

Update: I just added a little unit test that calls ElastixMain::GetComponentDatabase() and ComponentDatabase::GetCreator inside an OpenMP parallel for-loop. If that test passed successfully, I think I'll just merge this pull request.