Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Remove OpenMP based code, ELASTIX_USE_OPENMP option and OpenMP compiler flags #1188

Merged
merged 3 commits into from
Jul 8, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
17 changes: 0 additions & 17 deletions CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -128,23 +128,6 @@ if(ELASTIX_USE_EIGEN)
add_definitions(-DELASTIX_USE_EIGEN)
endif()

#---------------------------------------------------------------------
# Find OpenMP
mark_as_advanced(ELASTIX_USE_OPENMP)
option(ELASTIX_USE_OPENMP "Use OpenMP to speed up certain computations." ON)

if(ELASTIX_USE_OPENMP)
find_package(OpenMP QUIET)
if(OPENMP_FOUND)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
# set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS}
# ${OpenMP_EXE_LINKER_FLAGS}")
include_directories(${OpenMP_CXX_INCLUDE_DIRS})
link_libraries(${OpenMP_CXX_LIBRARIES})
add_definitions(-DELASTIX_USE_OPENMP)
endif()
endif()

#---------------------------------------------------------------------
# Set single (build-tree) output directories for all executables and libraries.
# This makes it easier to create an elxUseFile.cmake, that users can
Expand Down
5 changes: 0 additions & 5 deletions Common/CostFunctions/itkAdvancedImageToImageMetric.h
Original file line number Diff line number Diff line change
Expand Up @@ -298,10 +298,6 @@ class ITK_TEMPLATE_EXPORT AdvancedImageToImageMetric : public ImageToImageMetric
void
Initialize() override;

/** Set number of threads to use for computations. */
void
SetNumberOfWorkUnits(ThreadIdType numberOfThreads);

/** Switch the function BeforeThreadedGetValueAndDerivative on or off. */
itkSetMacro(UseMetricSingleThreaded, bool);
itkGetConstReferenceMacro(UseMetricSingleThreaded, bool);
Expand Down Expand Up @@ -423,7 +419,6 @@ class ITK_TEMPLATE_EXPORT AdvancedImageToImageMetric : public ImageToImageMetric
/** Variables for multi-threading. */
bool m_UseMetricSingleThreaded{ true };
bool m_UseMultiThread{ false };
bool m_UseOpenMP{};

/** Helper structs that multi-threads the computation of
* the metric derivative using ITK threads.
Expand Down
33 changes: 0 additions & 33 deletions Common/CostFunctions/itkAdvancedImageToImageMetric.hxx
Original file line number Diff line number Diff line change
Expand Up @@ -24,10 +24,6 @@
#include "itkAdvancedRayCastInterpolateImageFunction.h"
#include "itkComputeImageExtremaFilter.h"

#ifdef ELASTIX_USE_OPENMP
# include <omp.h>
#endif

#include <algorithm> // For min.
#include <cassert>

Expand All @@ -49,41 +45,12 @@ AdvancedImageToImageMetric<TFixedImage, TMovingImage>::AdvancedImageToImageMetri
*/
this->SetComputeGradient(false);

/** OpenMP related. Switch to on when available */
#ifdef ELASTIX_USE_OPENMP
m_UseOpenMP = true;

const int nthreads = static_cast<int>(Self::GetNumberOfWorkUnits());
omp_set_num_threads(nthreads);
#else
m_UseOpenMP = false;
#endif

/** Initialize the m_ThreaderMetricParameters. */
m_ThreaderMetricParameters.st_Metric = this;

} // end Constructor


/**
* ********************* SetNumberOfWorkUnits ****************************
*/

template <class TFixedImage, class TMovingImage>
void
AdvancedImageToImageMetric<TFixedImage, TMovingImage>::SetNumberOfWorkUnits(ThreadIdType numberOfThreads)
{
// Note: This is a workaround for ITK5, which renamed NumberOfThreads
// to NumberOfWorkUnits
Superclass::SetNumberOfWorkUnits(numberOfThreads);

#ifdef ELASTIX_USE_OPENMP
const int nthreads = static_cast<int>(Self::GetNumberOfWorkUnits());
omp_set_num_threads(nthreads);
#endif
} // end SetNumberOfWorkUnits()


/**
* ********************* Initialize ****************************
*/
Expand Down
2 changes: 1 addition & 1 deletion Common/GTesting/elxElastixMainGTest.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -25,7 +25,7 @@
// Tests retrieving the component data base and a component creator in parallel.
GTEST_TEST(ElastixMain, GetComponentDatabaseAndCreatorInParallel)
{
#pragma omp parallel for
// TODO Make iterations like this run in parallel, for the sake of the test.
for (auto i = 0; i <= 9; ++i)
{
const auto creator = elx::ElastixMain::GetComponentDatabase().GetCreator("Elastix", 1);
Expand Down
1 change: 0 additions & 1 deletion Common/GTesting/itkAdvancedImageToImageMetricGTest.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -66,7 +66,6 @@ class TestMetric : public AdvancedImageToImageMetric<TImage, TImage>
EXPECT_EQ(AdvancedImageToImageMetricType::m_MovingImageMaxLimit, 1);
EXPECT_EQ(AdvancedImageToImageMetricType::m_UseMetricSingleThreaded, true);
EXPECT_EQ(AdvancedImageToImageMetricType::m_UseMultiThread, false);
EXPECT_EQ(AdvancedImageToImageMetricType::m_UseOpenMP, bool{ AdvancedImageToImageMetricType::m_UseOpenMP });
EXPECT_EQ(AdvancedImageToImageMetricType::m_ThreaderMetricParameters.st_Metric, this);
EXPECT_EQ(AdvancedImageToImageMetricType::m_ThreaderMetricParameters.st_DerivativePointer, nullptr);
EXPECT_EQ(AdvancedImageToImageMetricType::m_ThreaderMetricParameters.st_NormalizationFactor, 0.0);
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -555,42 +555,12 @@ ParzenWindowMutualInformationImageToImageMetric<TFixedImage, TMovingImage>::Afte
const ThreadIdType numberOfThreads = Self::GetNumberOfWorkUnits();

/** Accumulate derivatives. */
// compute single-threadedly
if (!Superclass::m_UseMultiThread && false) // force multi-threaded
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This if (!x && false) was introduced with commit 2fc2c49, "HUGE merge of the performance branch", November 14, 2013

{
derivative = Superclass::m_GetValueAndDerivativePerThreadVariables[0].st_Derivative;
for (ThreadIdType i = 1; i < numberOfThreads; ++i)
{
derivative += Superclass::m_GetValueAndDerivativePerThreadVariables[i].st_Derivative;
}
}
#ifdef ELASTIX_USE_OPENMP
// compute multi-threadedly with openmp
else if (false) // Superclass::m_UseOpenMP )
{
const int spaceDimension = static_cast<int>(this->GetNumberOfParameters());

# pragma omp parallel for
for (int j = 0; j < spaceDimension; ++j)
{
DerivativeValueType sum = Superclass::m_GetValueAndDerivativePerThreadVariables[0].st_Derivative[j];
for (ThreadIdType i = 1; i < numberOfThreads; ++i)
{
sum += Superclass::m_GetValueAndDerivativePerThreadVariables[i].st_Derivative[j];
}
derivative[j] = sum;
}
}
#endif
// compute multi-threadedly with itk threads
else
{
Superclass::m_ThreaderMetricParameters.st_DerivativePointer = derivative.begin();
Superclass::m_ThreaderMetricParameters.st_NormalizationFactor = 1.0;

this->m_Threader->SetSingleMethodAndExecute(this->AccumulateDerivativesThreaderCallback,
&(Superclass::m_ThreaderMetricParameters));
}

} // end AfterThreadedComputeDerivativeLowMemory()

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -61,13 +61,6 @@ AdvancedMeanSquaresMetric<TElastix>::BeforeEachResolution()
configuration.ReadParameter(useNormalization, "UseNormalization", BaseComponent::GetComponentLabel(), level, 0);
this->SetUseNormalization(useNormalization);

/** Select the use of an OpenMP implementation for GetValueAndDerivative. */
std::string useOpenMP = configuration.GetCommandLineArgument("-useOpenMP_SSD");
if (useOpenMP == "true")
{
this->SetUseOpenMP(true);
}

} // end BeforeEachResolution()


Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -154,13 +154,6 @@ class ITK_TEMPLATE_EXPORT AdvancedMeanSquaresImageToImageMetric
itkSetMacro(UseNormalization, bool);
itkGetConstMacro(UseNormalization, bool);

/** If the compiler supports OpenMP, this flag specifies whether
* or not to use it. For this metric we have an OpenMP variant for
* GetValueAndDerivative(). It is also used at other places.
* Note that MS Visual Studio and gcc support OpenMP.
*/
itkSetMacro(UseOpenMP, bool);

protected:
AdvancedMeanSquaresImageToImageMetric();
~AdvancedMeanSquaresImageToImageMetric() override = default;
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -23,10 +23,6 @@
#include "itkMersenneTwisterRandomVariateGenerator.h"
#include "itkComputeImageExtremaFilter.h"

#ifdef ELASTIX_USE_OPENMP
# include <omp.h>
#endif

namespace itk
{

Expand Down Expand Up @@ -671,42 +667,12 @@ AdvancedMeanSquaresImageToImageMetric<TFixedImage, TMovingImage>::AfterThreadedG
value *= normal_sum;

/** Accumulate derivatives. */
// compute single-threadedly
if (!Superclass::m_UseMultiThread && false) // force multi-threaded
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This if (!x && false) was also introduced with commit 2fc2c49, "HUGE merge of the performance branch", November 14, 2013

{
derivative = Superclass::m_GetValueAndDerivativePerThreadVariables[0].st_Derivative * normal_sum;
for (ThreadIdType i = 1; i < numberOfThreads; ++i)
{
derivative += Superclass::m_GetValueAndDerivativePerThreadVariables[i].st_Derivative * normal_sum;
}
}
// compute multi-threadedly with itk threads
else if (true) // force ITK threads !Superclass::m_UseOpenMP )
{
Superclass::m_ThreaderMetricParameters.st_DerivativePointer = derivative.begin();
Superclass::m_ThreaderMetricParameters.st_NormalizationFactor = 1.0 / normal_sum;

this->m_Threader->SetSingleMethodAndExecute(this->AccumulateDerivativesThreaderCallback,
&(Superclass::m_ThreaderMetricParameters));
}
#ifdef ELASTIX_USE_OPENMP
// compute multi-threadedly with openmp
else
{
const int spaceDimension = static_cast<int>(this->GetNumberOfParameters());

# pragma omp parallel for
for (int j = 0; j < spaceDimension; ++j)
{
DerivativeValueType sum{};
for (ThreadIdType i = 0; i < numberOfThreads; ++i)
{
sum += Superclass::m_GetValueAndDerivativePerThreadVariables[i].st_Derivative[j];
}
derivative[j] = sum * normal_sum;
}
}
#endif

} // end AfterThreadedGetValueAndDerivative()

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -20,10 +20,6 @@

#include "itkAdvancedNormalizedCorrelationImageToImageMetric.h"

#ifdef ELASTIX_USE_OPENMP
# include <omp.h>
#endif

#include <algorithm> // For min.
#include <cassert>

Expand Down Expand Up @@ -642,34 +638,8 @@ AdvancedNormalizedCorrelationImageToImageMetric<TFixedImage, TMovingImage>::Afte
value = sfm / denom;

/** Calculate the metric derivative. */
// single-threaded
if (!Superclass::m_UseMultiThread && false) // force multi-threaded
{
DerivativeType & derivativeF = this->m_CorrelationGetValueAndDerivativePerThreadVariables[0].st_DerivativeF;
DerivativeType & derivativeM = this->m_CorrelationGetValueAndDerivativePerThreadVariables[0].st_DerivativeM;
DerivativeType & differential = this->m_CorrelationGetValueAndDerivativePerThreadVariables[0].st_Differential;

for (ThreadIdType i = 1; i < numberOfThreads; ++i)
{
derivativeF += this->m_CorrelationGetValueAndDerivativePerThreadVariables[i].st_DerivativeF;
derivativeM += this->m_CorrelationGetValueAndDerivativePerThreadVariables[i].st_DerivativeM;
differential += this->m_CorrelationGetValueAndDerivativePerThreadVariables[i].st_Differential;
}
// force multi-threaded

const auto numberOfParameters = this->GetNumberOfParameters();

/** Subtract things from derivativeF and derivativeM. */
double diff, derF, derM;
for (unsigned int i = 0; i < numberOfParameters; ++i)
{
diff = differential[i];
derF = derivativeF[i] - (sf / N) * diff;
derM = derivativeM[i] - (sm / N) * diff;
derivative[i] = (derF - (sfm / smm) * derM) / denom;
}
}
else if (true) // force !Superclass::m_UseOpenMP ) // multi-threaded using ITK threads
{
MultiThreaderAccumulateDerivativeType userData;

userData.st_Metric = const_cast<Self *>(this);
Expand All @@ -680,38 +650,6 @@ AdvancedNormalizedCorrelationImageToImageMetric<TFixedImage, TMovingImage>::Afte
userData.st_DerivativePointer = derivative.begin();

this->m_Threader->SetSingleMethodAndExecute(AccumulateDerivativesThreaderCallback, &userData);
}
#ifdef ELASTIX_USE_OPENMP
// compute multi-threadedly with openmp
else
{
const int spaceDimension = static_cast<int>(this->GetNumberOfParameters());
const AccumulateType sf_N = sf / N;
const AccumulateType sm_N = sm / N;
const AccumulateType sfm_smm = sfm / smm;

# pragma omp parallel for
for (int j = 0; j < spaceDimension; ++j)
{
DerivativeValueType derivativeF = this->m_CorrelationGetValueAndDerivativePerThreadVariables[0].st_DerivativeF[j];
DerivativeValueType derivativeM = this->m_CorrelationGetValueAndDerivativePerThreadVariables[0].st_DerivativeM[j];
DerivativeValueType differential =
this->m_CorrelationGetValueAndDerivativePerThreadVariables[0].st_Differential[j];

for (ThreadIdType i = 1; i < numberOfThreads; ++i)
{
derivativeF += this->m_CorrelationGetValueAndDerivativePerThreadVariables[i].st_DerivativeF[j];
derivativeM += this->m_CorrelationGetValueAndDerivativePerThreadVariables[i].st_DerivativeM[j];
differential += this->m_CorrelationGetValueAndDerivativePerThreadVariables[i].st_Differential[j];
}

derivativeF -= sf_N * differential;
derivativeM -= sm_N * differential;

derivative[j] = (derivativeF - sfm_smm * derivativeM) / denom;
}
} // end OpenMP
#endif

} // end AfterThreadedGetValueAndDerivative()

Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -20,10 +20,6 @@

#include "itkTransformBendingEnergyPenaltyTerm.h"

#ifdef ELASTIX_USE_OPENMP
# include <omp.h>
#endif

namespace itk
{

Expand Down Expand Up @@ -570,7 +566,6 @@ TransformBendingEnergyPenaltyTerm<TFixedImage, TScalarType>::AfterThreadedGetVal

/** Accumulate derivatives. */
// it seems that multi-threaded adding is faster than single-threaded
// it seems that openmp is faster than itk threads
// compute single-threadedly
if (!Superclass::m_UseMultiThread)
{
Expand All @@ -582,7 +577,7 @@ TransformBendingEnergyPenaltyTerm<TFixedImage, TScalarType>::AfterThreadedGetVal
derivative /= static_cast<DerivativeValueType>(Superclass::m_NumberOfPixelsCounted);
}
// compute multi-threadedly with itk threads
else if (!Superclass::m_UseOpenMP || true) // force
else
{
Superclass::m_ThreaderMetricParameters.st_DerivativePointer = derivative.begin();
Superclass::m_ThreaderMetricParameters.st_NormalizationFactor =
Expand All @@ -591,26 +586,6 @@ TransformBendingEnergyPenaltyTerm<TFixedImage, TScalarType>::AfterThreadedGetVal
this->m_Threader->SetSingleMethodAndExecute(this->AccumulateDerivativesThreaderCallback,
&(Superclass::m_ThreaderMetricParameters));
}
#ifdef ELASTIX_USE_OPENMP
// compute multi-threadedly with openmp
else
{
const DerivativeValueType numPix = static_cast<DerivativeValueType>(Superclass::m_NumberOfPixelsCounted);
const int nthreads = static_cast<int>(numberOfThreads);
omp_set_num_threads(nthreads);
const int spaceDimension = static_cast<int>(this->GetNumberOfParameters());
# pragma omp parallel for
for (int j = 0; j < spaceDimension; ++j)
{
DerivativeValueType sum{};
for (ThreadIdType i = 0; i < numberOfThreads; ++i)
{
sum += Superclass::m_GetValueAndDerivativePerThreadVariables[i].st_Derivative[j];
}
derivative[j] = sum / numPix;
}
}
#endif

} // end AfterThreadedGetValueAndDerivative()

Expand Down
Loading
Loading