diff --git a/.clang-format b/.clang-format index 92c500268c08..89187cd67bd3 100644 --- a/.clang-format +++ b/.clang-format @@ -146,6 +146,16 @@ Standard: Cpp11 StatementMacros: - Q_UNUSED - QT_REQUIRE_VERSION + - GCC_PRAGMA_PUSH + - GCC_PRAGMA_POP + - GCC_SUPPRESS_Wfloat_equal + - GCC_SUPPRESS_Wformat_nonliteral + - CLANG_PRAGMA_PUSH + - CLANG_PRAGMA_POP + - CLANG_SUPPRESS_Wcpp14_extensions + - INTEL_PRAGMA_WARN_PUSH + - INTEL_PRAGMA_WARN_POP + - INTEL_SUPPRESS_warning_1292 TabWidth: 2 UseTab: Never ... diff --git a/.github/workflows/itk_dict.txt b/.github/workflows/itk_dict.txt index 5a69dae67705..40f96c0490b7 100644 --- a/.github/workflows/itk_dict.txt +++ b/.github/workflows/itk_dict.txt @@ -44,7 +44,6 @@ Crofton's Cuadra Cz DDataFrameobject -Dz Dasarathy Deref Dij @@ -55,6 +54,7 @@ Dxx Dyx Dyy Dyz +Dz Dzx Dzy Dzz @@ -74,13 +74,13 @@ Fletch Freesurfer Frenet GEAdw +GPUSuperclass's Gao Gibbsian +Gifti Golby Graphviz Greenleaf -Gifti -GPUSuperclass's Guillen Haar Hackathon @@ -95,11 +95,11 @@ Hilden Hiraki Hirschberg Hostname -Intelli IOregion ITKIOMeshGifti ITKImageToIplImage ITKOptimizersv +Intelli Irfan Irix Isocenter @@ -740,6 +740,7 @@ ponderation posn postcondition posteriori +pragmas preallocated precalculate precalculating diff --git a/Modules/Core/Common/include/itkConceptChecking.h b/Modules/Core/Common/include/itkConceptChecking.h index be4bb8583621..4f679e7753b6 100644 --- a/Modules/Core/Common/include/itkConceptChecking.h +++ b/Modules/Core/Common/include/itkConceptChecking.h @@ -304,10 +304,11 @@ struct EqualityComparable void constraints() { - CLANG_PRAGMA_PUSH - CLANG_SUPPRESS_Wfloat_equal Detail::RequireBooleanExpression(a == b); + GCC_PRAGMA_PUSH + GCC_SUPPRESS_Wfloat_equal + Detail::RequireBooleanExpression(a == b); Detail::RequireBooleanExpression(a != b); - CLANG_PRAGMA_POP + GCC_PRAGMA_POP } T1 a; @@ -331,10 +332,11 @@ struct Comparable Detail::RequireBooleanExpression(a > b); Detail::RequireBooleanExpression(a <= b); Detail::RequireBooleanExpression(a >= b); - CLANG_PRAGMA_PUSH - CLANG_SUPPRESS_Wfloat_equal Detail::RequireBooleanExpression(a == b); + GCC_PRAGMA_PUSH + GCC_SUPPRESS_Wfloat_equal + Detail::RequireBooleanExpression(a == b); Detail::RequireBooleanExpression(a != b); - CLANG_PRAGMA_POP + GCC_PRAGMA_POP } T1 a; diff --git a/Modules/Core/Common/include/itkMacro.h b/Modules/Core/Common/include/itkMacro.h index a01303407c29..df2c284540e2 100644 --- a/Modules/Core/Common/include/itkMacro.h +++ b/Modules/Core/Common/include/itkMacro.h @@ -99,18 +99,27 @@ namespace itk // to be quoted. #define ITK_PRAGMA(x) _Pragma(#x) -// The clang compiler has many useful non-default compiler warnings -// that tend to have a high false positive rate. -// The following set of defines allows us to suppress false positives -// and still track down suspicious code +// The GCC/Clang compilers have many useful non-default compiler warnings +// that tend to have a high false positive rate or are otherwise not always appropriate. +// The following set of defines allows us to suppress instances of said warnings. + +// For GCC and Clang (Clang also identifies itself as GCC, and supports these pragmas): +#if defined(__GNUC__) +# define GCC_PRAGMA_PUSH ITK_PRAGMA(GCC diagnostic push) +# define GCC_PRAGMA_POP ITK_PRAGMA(GCC diagnostic pop) +# define GCC_SUPPRESS_Wfloat_equal ITK_PRAGMA(GCC diagnostic ignored "-Wfloat-equal") +# define GCC_SUPPRESS_Wformat_nonliteral ITK_PRAGMA(GCC diagnostic ignored "-Wformat-nonliteral") +#else +# define GCC_PRAGMA_PUSH +# define GCC_PRAGMA_POP +# define GCC_SUPPRESS_Wfloat_equal +# define GCC_SUPPRESS_Wformat_nonliteral +#endif + +// For Clang only (and not GCC): #if defined(__clang__) && defined(__has_warning) # define CLANG_PRAGMA_PUSH ITK_PRAGMA(clang diagnostic push) # define CLANG_PRAGMA_POP ITK_PRAGMA(clang diagnostic pop) -# if __has_warning("-Wfloat-equal") -# define CLANG_SUPPRESS_Wfloat_equal ITK_PRAGMA(clang diagnostic ignored "-Wfloat-equal") -# else -# define CLANG_SUPPRESS_Wfloat_equal -# endif # if __has_warning("-Wc++14-extensions") # define CLANG_SUPPRESS_Wcpp14_extensions ITK_PRAGMA(clang diagnostic ignored "-Wc++14-extensions") # else @@ -119,7 +128,6 @@ namespace itk #else # define CLANG_PRAGMA_PUSH # define CLANG_PRAGMA_POP -# define CLANG_SUPPRESS_Wfloat_equal # define CLANG_SUPPRESS_Wcpp14_extensions #endif @@ -887,13 +895,13 @@ compilers. itkDebugMacro("setting input " #name " to " << _arg); \ const DecoratorType * oldInput = \ itkDynamicCastInDebugMode(this->ProcessObject::GetInput(#name)); \ - CLANG_PRAGMA_PUSH \ - CLANG_SUPPRESS_Wfloat_equal \ + GCC_PRAGMA_PUSH \ + GCC_SUPPRESS_Wfloat_equal \ if (oldInput && oldInput->Get() == _arg) \ { \ return; \ } \ - CLANG_PRAGMA_POP \ + GCC_PRAGMA_POP \ auto newInput = DecoratorType::New(); \ newInput->Set(_arg); \ this->Set##name##Input(newInput); \ @@ -991,17 +999,17 @@ compilers. /** Set built-in type or regular C++ type. Creates member Set"name"() (e.g., SetVisibility()); */ // clang-format off #define itkSetMacro(name, type) \ - virtual void Set##name(type _arg) \ + virtual void Set##name(type _arg) \ { \ itkDebugMacro("setting " #name " to " << _arg); \ - CLANG_PRAGMA_PUSH \ - CLANG_SUPPRESS_Wfloat_equal \ + GCC_PRAGMA_PUSH \ + GCC_SUPPRESS_Wfloat_equal \ if (this->m_##name != _arg) \ { \ - this->m_##name = std::move(_arg); \ + this->m_##name = std::move(_arg); \ this->Modified(); \ } \ - CLANG_PRAGMA_POP \ + GCC_PRAGMA_POP \ } \ ITK_MACROEND_NOOP_STATEMENT // clang-format on @@ -1089,14 +1097,14 @@ compilers. { \ const type temp_extrema = (_arg <= min ? min : (_arg >= max ? max : _arg)); \ itkDebugMacro("setting " << #name " to " << _arg); \ - CLANG_PRAGMA_PUSH \ - CLANG_SUPPRESS_Wfloat_equal \ + GCC_PRAGMA_PUSH \ + GCC_SUPPRESS_Wfloat_equal \ if (this->m_##name != temp_extrema) \ { \ this->m_##name = temp_extrema; \ this->Modified(); \ } \ - CLANG_PRAGMA_POP \ + GCC_PRAGMA_POP \ } \ ITK_MACROEND_NOOP_STATEMENT // clang-format on @@ -1217,13 +1225,13 @@ compilers. unsigned int i; \ for (i = 0; i < count; ++i) \ { \ - CLANG_PRAGMA_PUSH \ - CLANG_SUPPRESS_Wfloat_equal \ + GCC_PRAGMA_PUSH \ + GCC_SUPPRESS_Wfloat_equal \ if (data[i] != this->m_##name[i]) \ { \ break; \ } \ - CLANG_PRAGMA_POP \ + GCC_PRAGMA_POP \ } \ if (i < count) \ { \ diff --git a/Modules/Core/Common/include/itkMath.h b/Modules/Core/Common/include/itkMath.h index e1b47c5a6449..8e9fe85592cc 100644 --- a/Modules/Core/Common/include/itkMath.h +++ b/Modules/Core/Common/include/itkMath.h @@ -724,9 +724,10 @@ template inline bool ExactlyEquals(const TInput1 & x1, const TInput2 & x2) { - CLANG_PRAGMA_PUSH - CLANG_SUPPRESS_Wfloat_equal return x1 == x2; - CLANG_PRAGMA_POP + GCC_PRAGMA_PUSH + GCC_SUPPRESS_Wfloat_equal + return x1 == x2; + GCC_PRAGMA_POP } // The NotExactlyEquals function diff --git a/Modules/Core/Common/include/itkMathDetail.h b/Modules/Core/Common/include/itkMathDetail.h index c819e86ad95e..690e4a3cb8bb 100644 --- a/Modules/Core/Common/include/itkMathDetail.h +++ b/Modules/Core/Common/include/itkMathDetail.h @@ -86,12 +86,12 @@ namespace Detail //////////////////////////////////////// // Base versions -CLANG_PRAGMA_PUSH -CLANG_SUPPRESS_Wfloat_equal +GCC_PRAGMA_PUSH +GCC_SUPPRESS_Wfloat_equal - template - inline TReturn - RoundHalfIntegerToEven_base(TInput x) +template +inline TReturn +RoundHalfIntegerToEven_base(TInput x) { if (NumericTraits::IsNonnegative(x)) { @@ -135,7 +135,7 @@ Ceil_base(TInput x) return (NumericTraits::IsNegative(x)) ? r : (x == static_cast(r) ? r : r + static_cast(1)); } -CLANG_PRAGMA_POP +GCC_PRAGMA_POP //////////////////////////////////////// // 32 bits versions diff --git a/Modules/Core/Common/include/itkMultiThreaderBase.h b/Modules/Core/Common/include/itkMultiThreaderBase.h index 5f54d646cbe1..81cca72f39ac 100644 --- a/Modules/Core/Common/include/itkMultiThreaderBase.h +++ b/Modules/Core/Common/include/itkMultiThreaderBase.h @@ -242,11 +242,11 @@ CLANG_PRAGMA_PUSH CLANG_SUPPRESS_Wcpp14_extensions // clang-format on # ifdef ITK_LEGACY_SILENT - struct ThreadInfoStruct + struct ThreadInfoStruct # else - struct [[deprecated("Use WorkUnitInfo, ThreadInfoStruct is deprecated since ITK 5.0")]] ThreadInfoStruct + struct [[deprecated("Use WorkUnitInfo, ThreadInfoStruct is deprecated since ITK 5.0")]] ThreadInfoStruct # endif - // clang-format off + // clang-format off CLANG_PRAGMA_POP INTEL_PRAGMA_WARN_POP // clang-format on diff --git a/Modules/Core/Common/test/itkMathRoundProfileTest1.cxx b/Modules/Core/Common/test/itkMathRoundProfileTest1.cxx index 88ef05e93562..ea51d4752efa 100644 --- a/Modules/Core/Common/test/itkMathRoundProfileTest1.cxx +++ b/Modules/Core/Common/test/itkMathRoundProfileTest1.cxx @@ -26,20 +26,25 @@ itkMathRoundTestHelperFunction(double x) return static_cast(x >= 0. ? x : (itk::Math::ExactlyEquals(x, static_cast(x)) ? x : x - 1.)); } -#define itkRoundMacro(x, y) \ - if (x >= 0.5) \ - { \ - y = static_cast(x + 0.5); \ - } \ - else \ - { \ - CLANG_PRAGMA_PUSH \ - CLANG_SUPPRESS_Wfloat_equal if ((x + 0.5) == static_cast(x + 0.5)) CLANG_PRAGMA_POP \ - { \ - y = static_cast(x + 0.5); \ - } \ - else { y = static_cast(x - 0.5); } \ - } \ +#define itkRoundMacro(x, y) \ + if (x >= 0.5) \ + { \ + y = static_cast(x + 0.5); \ + } \ + else \ + { \ + GCC_PRAGMA_PUSH \ + GCC_SUPPRESS_Wfloat_equal \ + if ((x + 0.5) == static_cast(x + 0.5)) \ + { \ + y = static_cast(x + 0.5); \ + } \ + else \ + { \ + y = static_cast(x - 0.5); \ + } \ + GCC_PRAGMA_POP \ + } \ ITK_MACROEND_NOOP_STATEMENT int diff --git a/Modules/Core/Common/test/itkRealTimeIntervalTest.cxx b/Modules/Core/Common/test/itkRealTimeIntervalTest.cxx index c0bd620820ad..4551371d6195 100644 --- a/Modules/Core/Common/test/itkRealTimeIntervalTest.cxx +++ b/Modules/Core/Common/test/itkRealTimeIntervalTest.cxx @@ -25,9 +25,10 @@ #define CHECK_FOR_VALUE(a, b) \ { \ double eps = 4.0 * itk::NumericTraits::epsilon(); \ - CLANG_PRAGMA_PUSH \ - CLANG_SUPPRESS_Wfloat_equal eps = (b == 0.0) ? eps : itk::Math::abs(b * eps); \ - CLANG_PRAGMA_POP \ + GCC_PRAGMA_PUSH \ + GCC_SUPPRESS_Wfloat_equal \ + eps = (b == 0.0) ? eps : itk::Math::abs(b * eps); \ + GCC_PRAGMA_POP \ if (itk::Math::abs(a - b) > eps) \ { \ std::cerr << "Error in " #a << " expected " << b << " but got " << a << std::endl; \ diff --git a/Modules/Core/Common/test/itkRealTimeStampTest.cxx b/Modules/Core/Common/test/itkRealTimeStampTest.cxx index 558213e049da..e84a3febe6ff 100644 --- a/Modules/Core/Common/test/itkRealTimeStampTest.cxx +++ b/Modules/Core/Common/test/itkRealTimeStampTest.cxx @@ -25,9 +25,10 @@ #define CHECK_FOR_VALUE(a, b) \ { \ double eps = 4.0 * itk::NumericTraits::epsilon(); \ - CLANG_PRAGMA_PUSH \ - CLANG_SUPPRESS_Wfloat_equal eps = (b == 0.0) ? eps : itk::Math::abs(b * eps); \ - CLANG_PRAGMA_POP \ + GCC_PRAGMA_PUSH \ + GCC_SUPPRESS_Wfloat_equal \ + eps = (b == 0.0) ? eps : itk::Math::abs(b * eps); \ + GCC_PRAGMA_POP \ if (itk::Math::abs(a - b) > eps) \ { \ std::cerr << "Error in " #a << " expected " << b << " but got " << a << std::endl; \ diff --git a/Modules/Core/Common/test/itkVariableLengthVectorTest.cxx b/Modules/Core/Common/test/itkVariableLengthVectorTest.cxx index b14cf04af13a..e31d6f59772d 100644 --- a/Modules/Core/Common/test/itkVariableLengthVectorTest.cxx +++ b/Modules/Core/Common/test/itkVariableLengthVectorTest.cxx @@ -21,12 +21,14 @@ #include "itkMath.h" #define ASSERT(cond, text) \ - CLANG_PRAGMA_PUSH \ - CLANG_SUPPRESS_Wfloat_equal if (!(cond)) CLANG_PRAGMA_POP \ + GCC_PRAGMA_PUSH \ + GCC_SUPPRESS_Wfloat_equal \ + if (!(cond)) \ { \ std::cerr << __FILE__ << ':' << __LINE__ << ':' << "Assertion failed: " << #cond << ": " << text << std::endl; \ result = EXIT_FAILURE; \ } \ + GCC_PRAGMA_POP \ ITK_MACROEND_NOOP_STATEMENT int diff --git a/Modules/Core/TestKernel/include/itkTestingMacros.h b/Modules/Core/TestKernel/include/itkTestingMacros.h index 4dd87c96da7b..07b374f51381 100644 --- a/Modules/Core/TestKernel/include/itkTestingMacros.h +++ b/Modules/Core/TestKernel/include/itkTestingMacros.h @@ -143,9 +143,10 @@ namespace itk #define ITK_TEST_EXPECT_TRUE_STATUS_VALUE(command, statusVal) \ { \ - CLANG_PRAGMA_PUSH \ - CLANG_SUPPRESS_Wfloat_equal bool _ITK_TEST_EXPECT_TRUE_command(command); \ - CLANG_PRAGMA_POP \ + GCC_PRAGMA_PUSH \ + GCC_SUPPRESS_Wfloat_equal \ + bool _ITK_TEST_EXPECT_TRUE_command(command); \ + GCC_PRAGMA_POP \ if (!(_ITK_TEST_EXPECT_TRUE_command)) \ { \ std::cerr << "Error in " << #command << std::endl; \ @@ -159,9 +160,10 @@ namespace itk #define ITK_TEST_EXPECT_TRUE(command) \ { \ - CLANG_PRAGMA_PUSH \ - CLANG_SUPPRESS_Wfloat_equal bool _ITK_TEST_EXPECT_TRUE_command(command); \ - CLANG_PRAGMA_POP \ + GCC_PRAGMA_PUSH \ + GCC_SUPPRESS_Wfloat_equal \ + bool _ITK_TEST_EXPECT_TRUE_command(command); \ + GCC_PRAGMA_POP \ if (!(_ITK_TEST_EXPECT_TRUE_command)) \ { \ std::cerr << "Error in " << #command << std::endl; \ @@ -174,38 +176,40 @@ namespace itk ITK_MACROEND_NOOP_STATEMENT -#define ITK_TEST_EXPECT_EQUAL_STATUS_VALUE(lh, rh, statusVal) \ - { \ - CLANG_PRAGMA_PUSH \ - CLANG_SUPPRESS_Wfloat_equal bool _ITK_TEST_EXPECT_EQUAL_result((lh) == (rh)); \ - CLANG_PRAGMA_POP \ - if (!(_ITK_TEST_EXPECT_EQUAL_result)) \ - { \ - std::cerr << "Error in " << #lh << " == " << #rh << std::endl; \ - std::cerr << "\tIn " __FILE__ ", line " << __LINE__ << std::endl; \ - std::cerr << "\tlh: " << (lh) << std::endl; \ - std::cerr << "\trh: " << (rh) << std::endl; \ - std::cerr << "Expression is not equal" << std::endl; \ - statusVal = EXIT_FAILURE; \ - } \ - } \ +#define ITK_TEST_EXPECT_EQUAL_STATUS_VALUE(lh, rh, statusVal) \ + { \ + GCC_PRAGMA_PUSH \ + GCC_SUPPRESS_Wfloat_equal \ + bool _ITK_TEST_EXPECT_EQUAL_result((lh) == (rh)); \ + GCC_PRAGMA_POP \ + if (!(_ITK_TEST_EXPECT_EQUAL_result)) \ + { \ + std::cerr << "Error in " << #lh << " == " << #rh << std::endl; \ + std::cerr << "\tIn " __FILE__ ", line " << __LINE__ << std::endl; \ + std::cerr << "\tlh: " << (lh) << std::endl; \ + std::cerr << "\trh: " << (rh) << std::endl; \ + std::cerr << "Expression is not equal" << std::endl; \ + statusVal = EXIT_FAILURE; \ + } \ + } \ ITK_MACROEND_NOOP_STATEMENT -#define ITK_TEST_EXPECT_EQUAL(lh, rh) \ - { \ - CLANG_PRAGMA_PUSH \ - CLANG_SUPPRESS_Wfloat_equal bool _ITK_TEST_EXPECT_EQUAL_result((lh) == (rh)); \ - CLANG_PRAGMA_POP \ - if (!(_ITK_TEST_EXPECT_EQUAL_result)) \ - { \ - std::cerr << "Error in " << #lh << " == " << #rh << std::endl; \ - std::cerr << "\tIn " __FILE__ ", line " << __LINE__ << std::endl; \ - std::cerr << "\tlh: " << (lh) << std::endl; \ - std::cerr << "\trh: " << (rh) << std::endl; \ - std::cerr << "Expression is not equal" << std::endl; \ - return EXIT_FAILURE; \ - } \ - } \ +#define ITK_TEST_EXPECT_EQUAL(lh, rh) \ + { \ + GCC_PRAGMA_PUSH \ + GCC_SUPPRESS_Wfloat_equal \ + bool _ITK_TEST_EXPECT_EQUAL_result((lh) == (rh)); \ + GCC_PRAGMA_POP \ + if (!(_ITK_TEST_EXPECT_EQUAL_result)) \ + { \ + std::cerr << "Error in " << #lh << " == " << #rh << std::endl; \ + std::cerr << "\tIn " __FILE__ ", line " << __LINE__ << std::endl; \ + std::cerr << "\tlh: " << (lh) << std::endl; \ + std::cerr << "\trh: " << (rh) << std::endl; \ + std::cerr << "Expression is not equal" << std::endl; \ + return EXIT_FAILURE; \ + } \ + } \ ITK_MACROEND_NOOP_STATEMENT @@ -221,16 +225,18 @@ namespace itk ITK_MACROEND_NOOP_STATEMENT -#define ITK_TEST_SET_GET_VALUE(variable, command) \ - CLANG_PRAGMA_PUSH \ - CLANG_SUPPRESS_Wfloat_equal if (variable != command) CLANG_PRAGMA_POP \ - { \ - std::cerr << "Error in " << #command << std::endl; \ - std::cerr << " In " __FILE__ ", line " << __LINE__ << std::endl; \ - std::cerr << "Expected " << variable << std::endl; \ - std::cerr << "but got " << command << std::endl; \ - return EXIT_FAILURE; \ - } \ +#define ITK_TEST_SET_GET_VALUE(variable, command) \ + GCC_PRAGMA_PUSH \ + GCC_SUPPRESS_Wfloat_equal \ + if (variable != command) \ + { \ + std::cerr << "Error in " << #command << std::endl; \ + std::cerr << " In " __FILE__ ", line " << __LINE__ << std::endl; \ + std::cerr << "Expected " << variable << std::endl; \ + std::cerr << "but got " << command << std::endl; \ + return EXIT_FAILURE; \ + } \ + GCC_PRAGMA_POP \ ITK_MACROEND_NOOP_STATEMENT #define ITK_TEST_SET_GET_NULL_VALUE(command) \ diff --git a/Modules/Filtering/BiasCorrection/include/itkN4BiasFieldCorrectionImageFilter.hxx b/Modules/Filtering/BiasCorrection/include/itkN4BiasFieldCorrectionImageFilter.hxx index cbb01ab28d13..8e50d1278b52 100644 --- a/Modules/Filtering/BiasCorrection/include/itkN4BiasFieldCorrectionImageFilter.hxx +++ b/Modules/Filtering/BiasCorrection/include/itkN4BiasFieldCorrectionImageFilter.hxx @@ -32,679 +32,684 @@ #include "itkSubtractImageFilter.h" #include "itkVectorIndexSelectionCastImageFilter.h" -CLANG_PRAGMA_PUSH -CLANG_SUPPRESS_Wfloat_equal +GCC_PRAGMA_PUSH +GCC_SUPPRESS_Wfloat_equal #include "vnl/algo/vnl_fft_1d.h" #include "vnl/vnl_complex_traits.h" #include "complex" - CLANG_PRAGMA_POP +GCC_PRAGMA_POP - namespace itk +namespace itk { - template - N4BiasFieldCorrectionImageFilter::N4BiasFieldCorrectionImageFilter() - : m_MaskLabel(NumericTraits::OneValue()) - , m_CurrentConvergenceMeasurement(RealType{}) +template +N4BiasFieldCorrectionImageFilter::N4BiasFieldCorrectionImageFilter() + : m_MaskLabel(NumericTraits::OneValue()) + , m_CurrentConvergenceMeasurement(RealType{}) - { - // implicit: - // #0 "Primary" required +{ + // implicit: + // #0 "Primary" required - // #1 "MaskImage" optional - Self::AddOptionalInputName("MaskImage", 1); + // #1 "MaskImage" optional + Self::AddOptionalInputName("MaskImage", 1); - // #2 "ConfidenceImage" optional - Self::AddOptionalInputName("ConfidenceImage", 2); + // #2 "ConfidenceImage" optional + Self::AddOptionalInputName("ConfidenceImage", 2); - this->SetNumberOfRequiredInputs(1); + this->SetNumberOfRequiredInputs(1); - this->m_LogBiasFieldControlPointLattice = nullptr; + this->m_LogBiasFieldControlPointLattice = nullptr; - this->m_NumberOfFittingLevels.Fill(1); - this->m_NumberOfControlPoints.Fill(4); + this->m_NumberOfFittingLevels.Fill(1); + this->m_NumberOfControlPoints.Fill(4); - this->m_MaximumNumberOfIterations.SetSize(1); - this->m_MaximumNumberOfIterations.Fill(50); - } + this->m_MaximumNumberOfIterations.SetSize(1); + this->m_MaximumNumberOfIterations.Fill(50); +} - template - void N4BiasFieldCorrectionImageFilter::EnlargeOutputRequestedRegion( - DataObject * output) - { - Superclass::EnlargeOutputRequestedRegion(output); +template +void +N4BiasFieldCorrectionImageFilter::EnlargeOutputRequestedRegion( + DataObject * output) +{ + Superclass::EnlargeOutputRequestedRegion(output); - if (output != nullptr) - { - output->SetRequestedRegionToLargestPossibleRegion(); - } + if (output != nullptr) + { + output->SetRequestedRegionToLargestPossibleRegion(); } +} - template - void N4BiasFieldCorrectionImageFilter::GenerateData() - { - this->AllocateOutputs(); +template +void +N4BiasFieldCorrectionImageFilter::GenerateData() +{ + this->AllocateOutputs(); - const InputImageType * inputImage = this->GetInput(); - using RegionType = typename InputImageType::RegionType; - const RegionType inputRegion = inputImage->GetBufferedRegion(); + const InputImageType * inputImage = this->GetInput(); + using RegionType = typename InputImageType::RegionType; + const RegionType inputRegion = inputImage->GetBufferedRegion(); - const typename InputImageType::SizeType inputImageSize = inputRegion.GetSize(); + const typename InputImageType::SizeType inputImageSize = inputRegion.GetSize(); - const MaskImageType * const maskImage = GetMaskImage(); + const MaskImageType * const maskImage = GetMaskImage(); - if ((maskImage != nullptr) && (maskImage->GetBufferedRegion().GetSize() != inputImageSize)) - { - itkExceptionMacro("If a mask image is specified, its size should be equal to the input image size"); - } + if ((maskImage != nullptr) && (maskImage->GetBufferedRegion().GetSize() != inputImageSize)) + { + itkExceptionMacro("If a mask image is specified, its size should be equal to the input image size"); + } - const RealImageType * const confidenceImage = GetConfidenceImage(); + const RealImageType * const confidenceImage = GetConfidenceImage(); - if ((confidenceImage != nullptr) && (confidenceImage->GetBufferedRegion().GetSize() != inputImageSize)) - { - itkExceptionMacro("If a confidence image is specified, its size should be equal to the input image size"); - } + if ((confidenceImage != nullptr) && (confidenceImage->GetBufferedRegion().GetSize() != inputImageSize)) + { + itkExceptionMacro("If a confidence image is specified, its size should be equal to the input image size"); + } - // Calculate the log of the input image. - RealImagePointer logInputImage = RealImageType::New(); - logInputImage->CopyInformation(inputImage); - logInputImage->SetRegions(inputRegion); - logInputImage->Allocate(false); + // Calculate the log of the input image. + RealImagePointer logInputImage = RealImageType::New(); + logInputImage->CopyInformation(inputImage); + logInputImage->SetRegions(inputRegion); + logInputImage->Allocate(false); - ImageAlgorithm::Copy(inputImage, logInputImage.GetPointer(), inputRegion, inputRegion); + ImageAlgorithm::Copy(inputImage, logInputImage.GetPointer(), inputRegion, inputRegion); - const auto maskImageBufferRange = MakeImageBufferRange(maskImage); - const auto confidenceImageBufferRange = MakeImageBufferRange(confidenceImage); - const MaskPixelType maskLabel = this->GetMaskLabel(); - const bool useMaskLabel = this->GetUseMaskLabel(); + const auto maskImageBufferRange = MakeImageBufferRange(maskImage); + const auto confidenceImageBufferRange = MakeImageBufferRange(confidenceImage); + const MaskPixelType maskLabel = this->GetMaskLabel(); + const bool useMaskLabel = this->GetUseMaskLabel(); - const ImageBufferRange logInputImageBufferRange{ *logInputImage }; - const size_t numberOfPixels = logInputImageBufferRange.size(); + const ImageBufferRange logInputImageBufferRange{ *logInputImage }; + const size_t numberOfPixels = logInputImageBufferRange.size(); - // Number of pixels of the input image that are included with the filter. - size_t numberOfIncludedPixels = 0; + // Number of pixels of the input image that are included with the filter. + size_t numberOfIncludedPixels = 0; - for (size_t indexValue = 0; indexValue < numberOfPixels; ++indexValue) + for (size_t indexValue = 0; indexValue < numberOfPixels; ++indexValue) + { + if ((maskImageBufferRange.empty() || (useMaskLabel && maskImageBufferRange[indexValue] == maskLabel) || + (!useMaskLabel && maskImageBufferRange[indexValue] != MaskPixelType{})) && + (confidenceImageBufferRange.empty() || confidenceImageBufferRange[indexValue] > 0.0)) { - if ((maskImageBufferRange.empty() || (useMaskLabel && maskImageBufferRange[indexValue] == maskLabel) || - (!useMaskLabel && maskImageBufferRange[indexValue] != MaskPixelType{})) && - (confidenceImageBufferRange.empty() || confidenceImageBufferRange[indexValue] > 0.0)) - { - ++numberOfIncludedPixels; - auto && logInputPixel = logInputImageBufferRange[indexValue]; + ++numberOfIncludedPixels; + auto && logInputPixel = logInputImageBufferRange[indexValue]; - if (logInputPixel > typename InputImageType::PixelType{}) - { - logInputPixel = std::log(static_cast(logInputPixel)); - } + if (logInputPixel > typename InputImageType::PixelType{}) + { + logInputPixel = std::log(static_cast(logInputPixel)); } } + } - // Duplicate logInputImage since we reuse the original at each iteration. + // Duplicate logInputImage since we reuse the original at each iteration. - using DuplicatorType = ImageDuplicator; - auto duplicator = DuplicatorType::New(); - duplicator->SetInputImage(logInputImage); - duplicator->Update(); + using DuplicatorType = ImageDuplicator; + auto duplicator = DuplicatorType::New(); + duplicator->SetInputImage(logInputImage); + duplicator->Update(); - RealImagePointer logUncorrectedImage = duplicator->GetOutput(); + RealImagePointer logUncorrectedImage = duplicator->GetOutput(); - // Provide an initial log bias field of zeros + // Provide an initial log bias field of zeros - RealImagePointer logBiasField = RealImageType::New(); - logBiasField->CopyInformation(inputImage); - logBiasField->SetRegions(inputImage->GetLargestPossibleRegion()); - logBiasField->AllocateInitialized(); + RealImagePointer logBiasField = RealImageType::New(); + logBiasField->CopyInformation(inputImage); + logBiasField->SetRegions(inputImage->GetLargestPossibleRegion()); + logBiasField->AllocateInitialized(); - RealImagePointer logSharpenedImage = RealImageType::New(); - logSharpenedImage->CopyInformation(inputImage); - logSharpenedImage->SetRegions(inputImage->GetLargestPossibleRegion()); - logSharpenedImage->Allocate(false); + RealImagePointer logSharpenedImage = RealImageType::New(); + logSharpenedImage->CopyInformation(inputImage); + logSharpenedImage->SetRegions(inputImage->GetLargestPossibleRegion()); + logSharpenedImage->Allocate(false); - // Iterate until convergence or iterative exhaustion. - unsigned int maximumNumberOfLevels = 1; - for (unsigned int d = 0; d < this->m_NumberOfFittingLevels.Size(); ++d) - { - if (this->m_NumberOfFittingLevels[d] > maximumNumberOfLevels) - { - maximumNumberOfLevels = this->m_NumberOfFittingLevels[d]; - } - } - if (this->m_MaximumNumberOfIterations.Size() != maximumNumberOfLevels) + // Iterate until convergence or iterative exhaustion. + unsigned int maximumNumberOfLevels = 1; + for (unsigned int d = 0; d < this->m_NumberOfFittingLevels.Size(); ++d) + { + if (this->m_NumberOfFittingLevels[d] > maximumNumberOfLevels) { - itkExceptionMacro("Number of iteration levels is not equal to the max number of levels."); + maximumNumberOfLevels = this->m_NumberOfFittingLevels[d]; } + } + if (this->m_MaximumNumberOfIterations.Size() != maximumNumberOfLevels) + { + itkExceptionMacro("Number of iteration levels is not equal to the max number of levels."); + } + + for (this->m_CurrentLevel = 0; this->m_CurrentLevel < maximumNumberOfLevels; this->m_CurrentLevel++) + { + IterationReporter reporter(this, 0, 1); - for (this->m_CurrentLevel = 0; this->m_CurrentLevel < maximumNumberOfLevels; this->m_CurrentLevel++) + this->m_ElapsedIterations = 0; + this->m_CurrentConvergenceMeasurement = NumericTraits::max(); + while (this->m_ElapsedIterations++ < this->m_MaximumNumberOfIterations[this->m_CurrentLevel] && + this->m_CurrentConvergenceMeasurement > this->m_ConvergenceThreshold) { - IterationReporter reporter(this, 0, 1); + // Sharpen the current estimate of the uncorrected image. + this->SharpenImage(logUncorrectedImage, logSharpenedImage); - this->m_ElapsedIterations = 0; - this->m_CurrentConvergenceMeasurement = NumericTraits::max(); - while (this->m_ElapsedIterations++ < this->m_MaximumNumberOfIterations[this->m_CurrentLevel] && - this->m_CurrentConvergenceMeasurement > this->m_ConvergenceThreshold) - { - // Sharpen the current estimate of the uncorrected image. - this->SharpenImage(logUncorrectedImage, logSharpenedImage); + using SubtracterType = SubtractImageFilter; + auto subtracter1 = SubtracterType::New(); + subtracter1->SetInput1(logUncorrectedImage); + subtracter1->SetInput2(logSharpenedImage); - using SubtracterType = SubtractImageFilter; - auto subtracter1 = SubtracterType::New(); - subtracter1->SetInput1(logUncorrectedImage); - subtracter1->SetInput2(logSharpenedImage); + RealImagePointer residualBiasField = subtracter1->GetOutput(); + residualBiasField->Update(); - RealImagePointer residualBiasField = subtracter1->GetOutput(); - residualBiasField->Update(); + // Smooth the residual bias field estimate and add the resulting + // control point grid to get the new total bias field estimate. - // Smooth the residual bias field estimate and add the resulting - // control point grid to get the new total bias field estimate. + RealImagePointer newLogBiasField = this->UpdateBiasFieldEstimate(residualBiasField, numberOfIncludedPixels); - RealImagePointer newLogBiasField = this->UpdateBiasFieldEstimate(residualBiasField, numberOfIncludedPixels); + this->m_CurrentConvergenceMeasurement = this->CalculateConvergenceMeasurement(logBiasField, newLogBiasField); + logBiasField = newLogBiasField; - this->m_CurrentConvergenceMeasurement = this->CalculateConvergenceMeasurement(logBiasField, newLogBiasField); - logBiasField = newLogBiasField; + auto subtracter2 = SubtracterType::New(); + subtracter2->SetInput1(logInputImage); + subtracter2->SetInput2(logBiasField); - auto subtracter2 = SubtracterType::New(); - subtracter2->SetInput1(logInputImage); - subtracter2->SetInput2(logBiasField); + logUncorrectedImage = subtracter2->GetOutput(); + logUncorrectedImage->Update(); - logUncorrectedImage = subtracter2->GetOutput(); - logUncorrectedImage->Update(); + reporter.CompletedStep(); + } - reporter.CompletedStep(); - } + using BSplineReconstructerType = BSplineControlPointImageFilter; + auto reconstructer = BSplineReconstructerType::New(); + reconstructer->SetInput(this->m_LogBiasFieldControlPointLattice); + reconstructer->SetOrigin(logBiasField->GetOrigin()); + reconstructer->SetSpacing(logBiasField->GetSpacing()); + reconstructer->SetDirection(logBiasField->GetDirection()); + reconstructer->SetSize(logBiasField->GetLargestPossibleRegion().GetSize()); + reconstructer->SetSplineOrder(this->m_SplineOrder); + reconstructer->Update(); - using BSplineReconstructerType = - BSplineControlPointImageFilter; - auto reconstructer = BSplineReconstructerType::New(); - reconstructer->SetInput(this->m_LogBiasFieldControlPointLattice); - reconstructer->SetOrigin(logBiasField->GetOrigin()); - reconstructer->SetSpacing(logBiasField->GetSpacing()); - reconstructer->SetDirection(logBiasField->GetDirection()); - reconstructer->SetSize(logBiasField->GetLargestPossibleRegion().GetSize()); - reconstructer->SetSplineOrder(this->m_SplineOrder); - reconstructer->Update(); - - typename BSplineReconstructerType::ArrayType numberOfLevels; - numberOfLevels.Fill(1); - for (unsigned int d = 0; d < ImageDimension; ++d) + typename BSplineReconstructerType::ArrayType numberOfLevels; + numberOfLevels.Fill(1); + for (unsigned int d = 0; d < ImageDimension; ++d) + { + if (this->m_NumberOfFittingLevels[d] + 1 >= this->m_CurrentLevel && + this->m_CurrentLevel != maximumNumberOfLevels - 1) { - if (this->m_NumberOfFittingLevels[d] + 1 >= this->m_CurrentLevel && - this->m_CurrentLevel != maximumNumberOfLevels - 1) - { - numberOfLevels[d] = 2; - } + numberOfLevels[d] = 2; } - this->m_LogBiasFieldControlPointLattice = reconstructer->RefineControlPointLattice(numberOfLevels); } - - using CustomBinaryFilter = itk::BinaryGeneratorImageFilter; - auto expAndDivFilter = CustomBinaryFilter::New(); - auto expAndDivLambda = - [](const typename InputImageType::PixelType & input, const typename RealImageType::PixelType & biasField) -> - typename OutputImageType::PixelType - { - return static_cast(input / std::exp(biasField)); - }; - expAndDivFilter->SetFunctor(expAndDivLambda); - expAndDivFilter->SetInput1(inputImage); - expAndDivFilter->SetInput2(logBiasField); - expAndDivFilter->Update(); - - this->GraftOutput(expAndDivFilter->GetOutput()); + this->m_LogBiasFieldControlPointLattice = reconstructer->RefineControlPointLattice(numberOfLevels); } - template - void N4BiasFieldCorrectionImageFilter::SharpenImage( - const RealImageType * unsharpenedImage, RealImageType * sharpenedImage) const + using CustomBinaryFilter = itk::BinaryGeneratorImageFilter; + auto expAndDivFilter = CustomBinaryFilter::New(); + auto expAndDivLambda = [](const typename InputImageType::PixelType & input, + const typename RealImageType::PixelType & biasField) -> typename OutputImageType::PixelType { - const auto maskImageBufferRange = MakeImageBufferRange(this->GetMaskImage()); - const auto confidenceImageBufferRange = MakeImageBufferRange(this->GetConfidenceImage()); - const MaskPixelType maskLabel = this->GetMaskLabel(); - const bool useMaskLabel = this->GetUseMaskLabel(); + return static_cast(input / std::exp(biasField)); + }; + expAndDivFilter->SetFunctor(expAndDivLambda); + expAndDivFilter->SetInput1(inputImage); + expAndDivFilter->SetInput2(logBiasField); + expAndDivFilter->Update(); + + this->GraftOutput(expAndDivFilter->GetOutput()); +} + +template +void +N4BiasFieldCorrectionImageFilter::SharpenImage( + const RealImageType * unsharpenedImage, + RealImageType * sharpenedImage) const +{ + const auto maskImageBufferRange = MakeImageBufferRange(this->GetMaskImage()); + const auto confidenceImageBufferRange = MakeImageBufferRange(this->GetConfidenceImage()); + const MaskPixelType maskLabel = this->GetMaskLabel(); + const bool useMaskLabel = this->GetUseMaskLabel(); - // Build the histogram for the uncorrected image. Store copy - // in a vnl_vector to utilize vnl FFT routines. Note that variables - // in real space are denoted by a single uppercase letter whereas their - // frequency counterparts are indicated by a trailing lowercase 'f'. + // Build the histogram for the uncorrected image. Store copy + // in a vnl_vector to utilize vnl FFT routines. Note that variables + // in real space are denoted by a single uppercase letter whereas their + // frequency counterparts are indicated by a trailing lowercase 'f'. - RealType binMaximum = NumericTraits::NonpositiveMin(); - RealType binMinimum = NumericTraits::max(); + RealType binMaximum = NumericTraits::NonpositiveMin(); + RealType binMinimum = NumericTraits::max(); - const auto unsharpenedImageBufferRange = MakeImageBufferRange(unsharpenedImage); - const size_t numberOfPixels = unsharpenedImageBufferRange.size(); + const auto unsharpenedImageBufferRange = MakeImageBufferRange(unsharpenedImage); + const size_t numberOfPixels = unsharpenedImageBufferRange.size(); - for (size_t indexValue = 0; indexValue < numberOfPixels; ++indexValue) + for (size_t indexValue = 0; indexValue < numberOfPixels; ++indexValue) + { + if ((maskImageBufferRange.empty() || (useMaskLabel && maskImageBufferRange[indexValue] == maskLabel) || + (!useMaskLabel && maskImageBufferRange[indexValue] != MaskPixelType{})) && + (confidenceImageBufferRange.empty() || confidenceImageBufferRange[indexValue] > 0.0)) { - if ((maskImageBufferRange.empty() || (useMaskLabel && maskImageBufferRange[indexValue] == maskLabel) || - (!useMaskLabel && maskImageBufferRange[indexValue] != MaskPixelType{})) && - (confidenceImageBufferRange.empty() || confidenceImageBufferRange[indexValue] > 0.0)) + RealType pixel = unsharpenedImageBufferRange[indexValue]; + if (pixel > binMaximum) + { + binMaximum = pixel; + } + else if (pixel < binMinimum) { - RealType pixel = unsharpenedImageBufferRange[indexValue]; - if (pixel > binMaximum) - { - binMaximum = pixel; - } - else if (pixel < binMinimum) - { - binMinimum = pixel; - } + binMinimum = pixel; } } - RealType histogramSlope = (binMaximum - binMinimum) / static_cast(this->m_NumberOfHistogramBins - 1); + } + RealType histogramSlope = (binMaximum - binMinimum) / static_cast(this->m_NumberOfHistogramBins - 1); - // Create the intensity profile (within the masked region, if applicable) - // using a triangular parzen windowing scheme. + // Create the intensity profile (within the masked region, if applicable) + // using a triangular parzen windowing scheme. - vnl_vector H(this->m_NumberOfHistogramBins, 0.0); + vnl_vector H(this->m_NumberOfHistogramBins, 0.0); - for (size_t indexValue = 0; indexValue < numberOfPixels; ++indexValue) + for (size_t indexValue = 0; indexValue < numberOfPixels; ++indexValue) + { + if ((maskImageBufferRange.empty() || (useMaskLabel && maskImageBufferRange[indexValue] == maskLabel) || + (!useMaskLabel && maskImageBufferRange[indexValue] != MaskPixelType{})) && + (confidenceImageBufferRange.empty() || confidenceImageBufferRange[indexValue] > 0.0)) { - if ((maskImageBufferRange.empty() || (useMaskLabel && maskImageBufferRange[indexValue] == maskLabel) || - (!useMaskLabel && maskImageBufferRange[indexValue] != MaskPixelType{})) && - (confidenceImageBufferRange.empty() || confidenceImageBufferRange[indexValue] > 0.0)) + RealType pixel = unsharpenedImageBufferRange[indexValue]; + + RealType cidx = (static_cast(pixel) - binMinimum) / histogramSlope; + unsigned int idx = itk::Math::floor(cidx); + RealType offset = cidx - static_cast(idx); + + if (offset == 0.0) { - RealType pixel = unsharpenedImageBufferRange[indexValue]; - - RealType cidx = (static_cast(pixel) - binMinimum) / histogramSlope; - unsigned int idx = itk::Math::floor(cidx); - RealType offset = cidx - static_cast(idx); - - if (offset == 0.0) - { - H[idx] += 1.0; - } - else if (idx < this->m_NumberOfHistogramBins - 1) - { - H[idx] += 1.0 - offset; - H[idx + 1] += offset; - } + H[idx] += 1.0; + } + else if (idx < this->m_NumberOfHistogramBins - 1) + { + H[idx] += 1.0 - offset; + H[idx + 1] += offset; } } + } - // Determine information about the intensity histogram and zero-pad - // histogram to a power of 2. + // Determine information about the intensity histogram and zero-pad + // histogram to a power of 2. - RealType exponent = std::ceil(std::log(static_cast(this->m_NumberOfHistogramBins)) / std::log(2.0)) + 1; - auto paddedHistogramSize = static_cast(std::pow(static_cast(2.0), exponent) + 0.5); - auto histogramOffset = static_cast(0.5 * (paddedHistogramSize - this->m_NumberOfHistogramBins)); + RealType exponent = std::ceil(std::log(static_cast(this->m_NumberOfHistogramBins)) / std::log(2.0)) + 1; + auto paddedHistogramSize = static_cast(std::pow(static_cast(2.0), exponent) + 0.5); + auto histogramOffset = static_cast(0.5 * (paddedHistogramSize - this->m_NumberOfHistogramBins)); - using FFTComputationType = double; - using FFTComplexType = std::complex; + using FFTComputationType = double; + using FFTComplexType = std::complex; - vnl_vector V(paddedHistogramSize, FFTComplexType(0.0, 0.0)); + vnl_vector V(paddedHistogramSize, FFTComplexType(0.0, 0.0)); - for (unsigned int n = 0; n < this->m_NumberOfHistogramBins; ++n) - { - V[n + histogramOffset] = H[n]; - } + for (unsigned int n = 0; n < this->m_NumberOfHistogramBins; ++n) + { + V[n + histogramOffset] = H[n]; + } - // Instantiate the 1-d vnl fft routine. + // Instantiate the 1-d vnl fft routine. - vnl_fft_1d fft(paddedHistogramSize); + vnl_fft_1d fft(paddedHistogramSize); - vnl_vector Vf(V); + vnl_vector Vf(V); - fft.fwd_transform(Vf); + fft.fwd_transform(Vf); - // Create the Gaussian filter. + // Create the Gaussian filter. - RealType scaledFWHM = this->m_BiasFieldFullWidthAtHalfMaximum / histogramSlope; - RealType expFactor = 4.0 * std::log(2.0) / itk::Math::sqr(scaledFWHM); - RealType scaleFactor = 2.0 * std::sqrt(std::log(2.0) / itk::Math::pi) / scaledFWHM; + RealType scaledFWHM = this->m_BiasFieldFullWidthAtHalfMaximum / histogramSlope; + RealType expFactor = 4.0 * std::log(2.0) / itk::Math::sqr(scaledFWHM); + RealType scaleFactor = 2.0 * std::sqrt(std::log(2.0) / itk::Math::pi) / scaledFWHM; - vnl_vector F(paddedHistogramSize, FFTComplexType(0.0, 0.0)); + vnl_vector F(paddedHistogramSize, FFTComplexType(0.0, 0.0)); - F[0] = FFTComplexType(scaleFactor, 0.0); - auto halfSize = static_cast(0.5 * paddedHistogramSize); - for (unsigned int n = 1; n <= halfSize; ++n) - { - F[n] = F[paddedHistogramSize - n] = - FFTComplexType(scaleFactor * std::exp(-itk::Math::sqr(static_cast(n)) * expFactor), 0.0); - } - if (paddedHistogramSize % 2 == 0) - { - F[halfSize] = FFTComplexType( - scaleFactor * std::exp(0.25 * -itk::Math::sqr(static_cast(paddedHistogramSize)) * expFactor), 0.0); - } + F[0] = FFTComplexType(scaleFactor, 0.0); + auto halfSize = static_cast(0.5 * paddedHistogramSize); + for (unsigned int n = 1; n <= halfSize; ++n) + { + F[n] = F[paddedHistogramSize - n] = + FFTComplexType(scaleFactor * std::exp(-itk::Math::sqr(static_cast(n)) * expFactor), 0.0); + } + if (paddedHistogramSize % 2 == 0) + { + F[halfSize] = FFTComplexType( + scaleFactor * std::exp(0.25 * -itk::Math::sqr(static_cast(paddedHistogramSize)) * expFactor), 0.0); + } - vnl_vector Ff(F); + vnl_vector Ff(F); - fft.fwd_transform(Ff); + fft.fwd_transform(Ff); - // Create the Wiener deconvolution filter. + // Create the Wiener deconvolution filter. - vnl_vector Gf(paddedHistogramSize); + vnl_vector Gf(paddedHistogramSize); - const auto wienerNoiseValue = static_cast(this->m_WienerFilterNoise); - for (unsigned int n = 0; n < paddedHistogramSize; ++n) - { - FFTComplexType c = vnl_complex_traits::conjugate(Ff[n]); - Gf[n] = c / (c * Ff[n] + wienerNoiseValue); - } + const auto wienerNoiseValue = static_cast(this->m_WienerFilterNoise); + for (unsigned int n = 0; n < paddedHistogramSize; ++n) + { + FFTComplexType c = vnl_complex_traits::conjugate(Ff[n]); + Gf[n] = c / (c * Ff[n] + wienerNoiseValue); + } - vnl_vector Uf(paddedHistogramSize); + vnl_vector Uf(paddedHistogramSize); - for (unsigned int n = 0; n < paddedHistogramSize; ++n) - { - Uf[n] = Vf[n] * Gf[n].real(); - } + for (unsigned int n = 0; n < paddedHistogramSize; ++n) + { + Uf[n] = Vf[n] * Gf[n].real(); + } - vnl_vector U(Uf); + vnl_vector U(Uf); - fft.bwd_transform(U); - for (unsigned int n = 0; n < paddedHistogramSize; ++n) - { - U[n] = FFTComplexType(std::max(U[n].real(), static_cast(0.0)), 0.0); - } + fft.bwd_transform(U); + for (unsigned int n = 0; n < paddedHistogramSize; ++n) + { + U[n] = FFTComplexType(std::max(U[n].real(), static_cast(0.0)), 0.0); + } - // Compute mapping E(u|v). + // Compute mapping E(u|v). - vnl_vector numerator(paddedHistogramSize); + vnl_vector numerator(paddedHistogramSize); - for (unsigned int n = 0; n < paddedHistogramSize; ++n) - { - numerator[n] = - FFTComplexType((binMinimum + (static_cast(n) - histogramOffset) * histogramSlope) * U[n].real(), 0.0); - } - fft.fwd_transform(numerator); - for (unsigned int n = 0; n < paddedHistogramSize; ++n) - { - numerator[n] *= Ff[n]; - } - fft.bwd_transform(numerator); + for (unsigned int n = 0; n < paddedHistogramSize; ++n) + { + numerator[n] = + FFTComplexType((binMinimum + (static_cast(n) - histogramOffset) * histogramSlope) * U[n].real(), 0.0); + } + fft.fwd_transform(numerator); + for (unsigned int n = 0; n < paddedHistogramSize; ++n) + { + numerator[n] *= Ff[n]; + } + fft.bwd_transform(numerator); - vnl_vector denominator(U); + vnl_vector denominator(U); - fft.fwd_transform(denominator); - for (unsigned int n = 0; n < paddedHistogramSize; ++n) - { - denominator[n] *= Ff[n]; - } - fft.bwd_transform(denominator); + fft.fwd_transform(denominator); + for (unsigned int n = 0; n < paddedHistogramSize; ++n) + { + denominator[n] *= Ff[n]; + } + fft.bwd_transform(denominator); - vnl_vector E(paddedHistogramSize); + vnl_vector E(paddedHistogramSize); - for (unsigned int n = 0; n < paddedHistogramSize; ++n) + for (unsigned int n = 0; n < paddedHistogramSize; ++n) + { + if (denominator[n].real() != 0.0) { - if (denominator[n].real() != 0.0) - { - E[n] = numerator[n].real() / denominator[n].real(); - } - else - { - E[n] = 0.0; - } + E[n] = numerator[n].real() / denominator[n].real(); } + else + { + E[n] = 0.0; + } + } - // Remove the zero-padding from the mapping. + // Remove the zero-padding from the mapping. - E = E.extract(this->m_NumberOfHistogramBins, histogramOffset); + E = E.extract(this->m_NumberOfHistogramBins, histogramOffset); - // Sharpen the image with the new mapping, E(u|v) - sharpenedImage->FillBuffer(0); + // Sharpen the image with the new mapping, E(u|v) + sharpenedImage->FillBuffer(0); - const ImageBufferRange sharpenedImageBufferRange{ *sharpenedImage }; + const ImageBufferRange sharpenedImageBufferRange{ *sharpenedImage }; - for (size_t indexValue = 0; indexValue < numberOfPixels; ++indexValue) + for (size_t indexValue = 0; indexValue < numberOfPixels; ++indexValue) + { + if ((maskImageBufferRange.empty() || (useMaskLabel && maskImageBufferRange[indexValue] == maskLabel) || + (!useMaskLabel && maskImageBufferRange[indexValue] != MaskPixelType{})) && + (confidenceImageBufferRange.empty() || confidenceImageBufferRange[indexValue] > 0.0)) { - if ((maskImageBufferRange.empty() || (useMaskLabel && maskImageBufferRange[indexValue] == maskLabel) || - (!useMaskLabel && maskImageBufferRange[indexValue] != MaskPixelType{})) && - (confidenceImageBufferRange.empty() || confidenceImageBufferRange[indexValue] > 0.0)) + RealType cidx = (unsharpenedImageBufferRange[indexValue] - binMinimum) / histogramSlope; + unsigned int idx = itk::Math::floor(cidx); + + RealType correctedPixel = 0; + if (idx < E.size() - 1) + { + correctedPixel = E[idx] + (E[idx + 1] - E[idx]) * (cidx - static_cast(idx)); + } + else { - RealType cidx = (unsharpenedImageBufferRange[indexValue] - binMinimum) / histogramSlope; - unsigned int idx = itk::Math::floor(cidx); - - RealType correctedPixel = 0; - if (idx < E.size() - 1) - { - correctedPixel = E[idx] + (E[idx + 1] - E[idx]) * (cidx - static_cast(idx)); - } - else - { - correctedPixel = E.back(); - } - sharpenedImageBufferRange[indexValue] = correctedPixel; + correctedPixel = E.back(); } + sharpenedImageBufferRange[indexValue] = correctedPixel; } } +} - template - typename N4BiasFieldCorrectionImageFilter::RealImagePointer - N4BiasFieldCorrectionImageFilter::UpdateBiasFieldEstimate( - RealImageType * fieldEstimate, const size_t numberOfIncludedPixels) +template +typename N4BiasFieldCorrectionImageFilter::RealImagePointer +N4BiasFieldCorrectionImageFilter::UpdateBiasFieldEstimate( + RealImageType * fieldEstimate, + const size_t numberOfIncludedPixels) +{ + // Temporarily set the direction cosine to identity since the B-spline + // approximation algorithm works in parametric space and not physical + // space. + typename ScalarImageType::DirectionType identity; + identity.SetIdentity(); + + const typename ScalarImageType::RegionType & bufferedRegion = fieldEstimate->GetBufferedRegion(); + const SizeValueType numberOfPixels = bufferedRegion.GetNumberOfPixels(); + const bool filterHandlesMemory = false; + + using ImporterType = ImportImageFilter; + auto importer = ImporterType::New(); + importer->SetImportPointer(fieldEstimate->GetBufferPointer(), numberOfPixels, filterHandlesMemory); + importer->SetRegion(fieldEstimate->GetBufferedRegion()); + importer->SetOrigin(fieldEstimate->GetOrigin()); + importer->SetSpacing(fieldEstimate->GetSpacing()); + importer->SetDirection(identity); + importer->Update(); + + const typename ImporterType::OutputImageType * parametricFieldEstimate = importer->GetOutput(); + + PointSetPointer fieldPoints = PointSetType::New(); + fieldPoints->Initialize(); + auto & pointSTLContainer = fieldPoints->GetPoints()->CastToSTLContainer(); + pointSTLContainer.reserve(numberOfIncludedPixels); + auto & pointDataSTLContainer = fieldPoints->GetPointData()->CastToSTLContainer(); + pointDataSTLContainer.reserve(numberOfIncludedPixels); + + typename BSplineFilterType::WeightsContainerType::Pointer weights = BSplineFilterType::WeightsContainerType::New(); + weights->Initialize(); + auto & weightSTLContainer = weights->CastToSTLContainer(); + weightSTLContainer.reserve(numberOfIncludedPixels); + + const auto maskImageBufferRange = MakeImageBufferRange(this->GetMaskImage()); + const auto confidenceImageBufferRange = MakeImageBufferRange(this->GetConfidenceImage()); + const MaskPixelType maskLabel = this->GetMaskLabel(); + const bool useMaskLabel = this->GetUseMaskLabel(); + + ImageRegionConstIteratorWithIndex It(parametricFieldEstimate, + parametricFieldEstimate->GetRequestedRegion()); + + for (size_t indexValue = 0; indexValue < numberOfPixels; ++indexValue, ++It) { - // Temporarily set the direction cosine to identity since the B-spline - // approximation algorithm works in parametric space and not physical - // space. - typename ScalarImageType::DirectionType identity; - identity.SetIdentity(); - - const typename ScalarImageType::RegionType & bufferedRegion = fieldEstimate->GetBufferedRegion(); - const SizeValueType numberOfPixels = bufferedRegion.GetNumberOfPixels(); - const bool filterHandlesMemory = false; - - using ImporterType = ImportImageFilter; - auto importer = ImporterType::New(); - importer->SetImportPointer(fieldEstimate->GetBufferPointer(), numberOfPixels, filterHandlesMemory); - importer->SetRegion(fieldEstimate->GetBufferedRegion()); - importer->SetOrigin(fieldEstimate->GetOrigin()); - importer->SetSpacing(fieldEstimate->GetSpacing()); - importer->SetDirection(identity); - importer->Update(); - - const typename ImporterType::OutputImageType * parametricFieldEstimate = importer->GetOutput(); - - PointSetPointer fieldPoints = PointSetType::New(); - fieldPoints->Initialize(); - auto & pointSTLContainer = fieldPoints->GetPoints()->CastToSTLContainer(); - pointSTLContainer.reserve(numberOfIncludedPixels); - auto & pointDataSTLContainer = fieldPoints->GetPointData()->CastToSTLContainer(); - pointDataSTLContainer.reserve(numberOfIncludedPixels); - - typename BSplineFilterType::WeightsContainerType::Pointer weights = BSplineFilterType::WeightsContainerType::New(); - weights->Initialize(); - auto & weightSTLContainer = weights->CastToSTLContainer(); - weightSTLContainer.reserve(numberOfIncludedPixels); - - const auto maskImageBufferRange = MakeImageBufferRange(this->GetMaskImage()); - const auto confidenceImageBufferRange = MakeImageBufferRange(this->GetConfidenceImage()); - const MaskPixelType maskLabel = this->GetMaskLabel(); - const bool useMaskLabel = this->GetUseMaskLabel(); - - ImageRegionConstIteratorWithIndex It(parametricFieldEstimate, - parametricFieldEstimate->GetRequestedRegion()); - - for (size_t indexValue = 0; indexValue < numberOfPixels; ++indexValue, ++It) + if ((maskImageBufferRange.empty() || (useMaskLabel && maskImageBufferRange[indexValue] == maskLabel) || + (!useMaskLabel && maskImageBufferRange[indexValue] != MaskPixelType{})) && + (confidenceImageBufferRange.empty() || confidenceImageBufferRange[indexValue] > 0.0)) { - if ((maskImageBufferRange.empty() || (useMaskLabel && maskImageBufferRange[indexValue] == maskLabel) || - (!useMaskLabel && maskImageBufferRange[indexValue] != MaskPixelType{})) && - (confidenceImageBufferRange.empty() || confidenceImageBufferRange[indexValue] > 0.0)) - { - PointType point; - parametricFieldEstimate->TransformIndexToPhysicalPoint(It.GetIndex(), point); + PointType point; + parametricFieldEstimate->TransformIndexToPhysicalPoint(It.GetIndex(), point); - ScalarType scalar; - scalar[0] = It.Get(); + ScalarType scalar; + scalar[0] = It.Get(); - pointDataSTLContainer.push_back(scalar); - pointSTLContainer.push_back(point); + pointDataSTLContainer.push_back(scalar); + pointSTLContainer.push_back(point); - RealType confidenceWeight = 1.0; - if (!confidenceImageBufferRange.empty()) - { - confidenceWeight = confidenceImageBufferRange[indexValue]; - } - weightSTLContainer.push_back(confidenceWeight); - } - } - - auto bspliner = BSplineFilterType::New(); - - typename BSplineFilterType::ArrayType numberOfControlPoints; - typename BSplineFilterType::ArrayType numberOfFittingLevels; - numberOfFittingLevels.Fill(1); - for (unsigned int d = 0; d < ImageDimension; ++d) - { - if (!this->m_LogBiasFieldControlPointLattice) + RealType confidenceWeight = 1.0; + if (!confidenceImageBufferRange.empty()) { - numberOfControlPoints[d] = this->m_NumberOfControlPoints[d]; - } - else - { - numberOfControlPoints[d] = this->m_LogBiasFieldControlPointLattice->GetLargestPossibleRegion().GetSize()[d]; + confidenceWeight = confidenceImageBufferRange[indexValue]; } + weightSTLContainer.push_back(confidenceWeight); } + } - typename ScalarImageType::PointType parametricOrigin = fieldEstimate->GetOrigin(); - for (unsigned int d = 0; d < ImageDimension; ++d) - { - parametricOrigin[d] += (fieldEstimate->GetSpacing()[d] * fieldEstimate->GetLargestPossibleRegion().GetIndex()[d]); - } - bspliner->SetOrigin(parametricOrigin); - bspliner->SetSpacing(fieldEstimate->GetSpacing()); - bspliner->SetSize(fieldEstimate->GetLargestPossibleRegion().GetSize()); - bspliner->SetDirection(fieldEstimate->GetDirection()); - bspliner->SetGenerateOutputImage(false); - bspliner->SetNumberOfLevels(numberOfFittingLevels); - bspliner->SetSplineOrder(this->m_SplineOrder); - bspliner->SetNumberOfControlPoints(numberOfControlPoints); - bspliner->SetInput(fieldPoints); - bspliner->SetPointWeights(weights); - bspliner->Update(); - - typename BiasFieldControlPointLatticeType::Pointer phiLattice = bspliner->GetPhiLattice(); - - // Add the bias field control points to the current estimate. + auto bspliner = BSplineFilterType::New(); + typename BSplineFilterType::ArrayType numberOfControlPoints; + typename BSplineFilterType::ArrayType numberOfFittingLevels; + numberOfFittingLevels.Fill(1); + for (unsigned int d = 0; d < ImageDimension; ++d) + { if (!this->m_LogBiasFieldControlPointLattice) { - this->m_LogBiasFieldControlPointLattice = phiLattice; + numberOfControlPoints[d] = this->m_NumberOfControlPoints[d]; } else { - // Ensure that the two lattices occupy the same physical space. Not - // necessary for performance since the parameters of the reconstructed - // bias field are specified later in this function in the reconstructer. - phiLattice->CopyInformation(this->m_LogBiasFieldControlPointLattice); - - using AdderType = AddImageFilter; - auto adder = AdderType::New(); - adder->SetInput1(this->m_LogBiasFieldControlPointLattice); - adder->SetInput2(phiLattice); - adder->Update(); - - this->m_LogBiasFieldControlPointLattice = adder->GetOutput(); + numberOfControlPoints[d] = this->m_LogBiasFieldControlPointLattice->GetLargestPossibleRegion().GetSize()[d]; } - - RealImagePointer smoothField = this->ReconstructBiasField(this->m_LogBiasFieldControlPointLattice); - smoothField->SetRegions(this->GetInput()->GetRequestedRegion()); - return smoothField; } - template - typename N4BiasFieldCorrectionImageFilter::RealImagePointer - N4BiasFieldCorrectionImageFilter::ReconstructBiasField( - const BiasFieldControlPointLatticeType * controlPointLattice) + typename ScalarImageType::PointType parametricOrigin = fieldEstimate->GetOrigin(); + for (unsigned int d = 0; d < ImageDimension; ++d) { - const InputImageType * inputImage = this->GetInput(); + parametricOrigin[d] += (fieldEstimate->GetSpacing()[d] * fieldEstimate->GetLargestPossibleRegion().GetIndex()[d]); + } + bspliner->SetOrigin(parametricOrigin); + bspliner->SetSpacing(fieldEstimate->GetSpacing()); + bspliner->SetSize(fieldEstimate->GetLargestPossibleRegion().GetSize()); + bspliner->SetDirection(fieldEstimate->GetDirection()); + bspliner->SetGenerateOutputImage(false); + bspliner->SetNumberOfLevels(numberOfFittingLevels); + bspliner->SetSplineOrder(this->m_SplineOrder); + bspliner->SetNumberOfControlPoints(numberOfControlPoints); + bspliner->SetInput(fieldPoints); + bspliner->SetPointWeights(weights); + bspliner->Update(); + + typename BiasFieldControlPointLatticeType::Pointer phiLattice = bspliner->GetPhiLattice(); + + // Add the bias field control points to the current estimate. + + if (!this->m_LogBiasFieldControlPointLattice) + { + this->m_LogBiasFieldControlPointLattice = phiLattice; + } + else + { + // Ensure that the two lattices occupy the same physical space. Not + // necessary for performance since the parameters of the reconstructed + // bias field are specified later in this function in the reconstructer. + phiLattice->CopyInformation(this->m_LogBiasFieldControlPointLattice); + + using AdderType = AddImageFilter; + auto adder = AdderType::New(); + adder->SetInput1(this->m_LogBiasFieldControlPointLattice); + adder->SetInput2(phiLattice); + adder->Update(); + + this->m_LogBiasFieldControlPointLattice = adder->GetOutput(); + } - using BSplineReconstructerType = BSplineControlPointImageFilter; - auto reconstructer = BSplineReconstructerType::New(); - reconstructer->SetInput(controlPointLattice); - reconstructer->SetOrigin(inputImage->GetOrigin()); - reconstructer->SetSpacing(inputImage->GetSpacing()); - reconstructer->SetDirection(inputImage->GetDirection()); - reconstructer->SetSplineOrder(this->m_SplineOrder); - reconstructer->SetSize(inputImage->GetLargestPossibleRegion().GetSize()); + RealImagePointer smoothField = this->ReconstructBiasField(this->m_LogBiasFieldControlPointLattice); + smoothField->SetRegions(this->GetInput()->GetRequestedRegion()); + return smoothField; +} - typename ScalarImageType::Pointer biasFieldBsplineImage = reconstructer->GetOutput(); - biasFieldBsplineImage->Update(); +template +typename N4BiasFieldCorrectionImageFilter::RealImagePointer +N4BiasFieldCorrectionImageFilter::ReconstructBiasField( + const BiasFieldControlPointLatticeType * controlPointLattice) +{ + const InputImageType * inputImage = this->GetInput(); + + using BSplineReconstructerType = BSplineControlPointImageFilter; + auto reconstructer = BSplineReconstructerType::New(); + reconstructer->SetInput(controlPointLattice); + reconstructer->SetOrigin(inputImage->GetOrigin()); + reconstructer->SetSpacing(inputImage->GetSpacing()); + reconstructer->SetDirection(inputImage->GetDirection()); + reconstructer->SetSplineOrder(this->m_SplineOrder); + reconstructer->SetSize(inputImage->GetLargestPossibleRegion().GetSize()); + + typename ScalarImageType::Pointer biasFieldBsplineImage = reconstructer->GetOutput(); + biasFieldBsplineImage->Update(); + + using SelectorType = VectorIndexSelectionCastImageFilter; + auto selector = SelectorType::New(); + selector->SetInput(biasFieldBsplineImage); + selector->SetIndex(0); + + RealImagePointer biasField = selector->GetOutput(); + biasField->Update(); + + biasField->DisconnectPipeline(); + + return biasField; +} + +template +typename N4BiasFieldCorrectionImageFilter::RealType +N4BiasFieldCorrectionImageFilter::CalculateConvergenceMeasurement( + const RealImageType * fieldEstimate1, + const RealImageType * fieldEstimate2) const +{ + using SubtracterType = SubtractImageFilter; + auto subtracter = SubtracterType::New(); + subtracter->SetInput1(fieldEstimate1); + subtracter->SetInput2(fieldEstimate2); + subtracter->Update(); - using SelectorType = VectorIndexSelectionCastImageFilter; - auto selector = SelectorType::New(); - selector->SetInput(biasFieldBsplineImage); - selector->SetIndex(0); + // Calculate statistics over the mask region - RealImagePointer biasField = selector->GetOutput(); - biasField->Update(); + RealType mu = 0.0; + RealType sigma = 0.0; + RealType N = 0.0; - biasField->DisconnectPipeline(); + const auto maskImageBufferRange = MakeImageBufferRange(this->GetMaskImage()); + const auto confidenceImageBufferRange = MakeImageBufferRange(this->GetConfidenceImage()); + const MaskPixelType maskLabel = this->GetMaskLabel(); + const bool useMaskLabel = this->GetUseMaskLabel(); - return biasField; - } + const auto subtracterImageBufferRange = MakeImageBufferRange(subtracter->GetOutput()); + const size_t numberOfPixels = subtracterImageBufferRange.size(); - template - typename N4BiasFieldCorrectionImageFilter::RealType - N4BiasFieldCorrectionImageFilter::CalculateConvergenceMeasurement( - const RealImageType * fieldEstimate1, const RealImageType * fieldEstimate2) const + for (size_t indexValue = 0; indexValue < numberOfPixels; ++indexValue) { - using SubtracterType = SubtractImageFilter; - auto subtracter = SubtracterType::New(); - subtracter->SetInput1(fieldEstimate1); - subtracter->SetInput2(fieldEstimate2); - subtracter->Update(); - - // Calculate statistics over the mask region - - RealType mu = 0.0; - RealType sigma = 0.0; - RealType N = 0.0; - - const auto maskImageBufferRange = MakeImageBufferRange(this->GetMaskImage()); - const auto confidenceImageBufferRange = MakeImageBufferRange(this->GetConfidenceImage()); - const MaskPixelType maskLabel = this->GetMaskLabel(); - const bool useMaskLabel = this->GetUseMaskLabel(); - - const auto subtracterImageBufferRange = MakeImageBufferRange(subtracter->GetOutput()); - const size_t numberOfPixels = subtracterImageBufferRange.size(); - - for (size_t indexValue = 0; indexValue < numberOfPixels; ++indexValue) + if ((maskImageBufferRange.empty() || (useMaskLabel && maskImageBufferRange[indexValue] == maskLabel) || + (!useMaskLabel && maskImageBufferRange[indexValue] != MaskPixelType{})) && + (confidenceImageBufferRange.empty() || confidenceImageBufferRange[indexValue] > 0.0)) { - if ((maskImageBufferRange.empty() || (useMaskLabel && maskImageBufferRange[indexValue] == maskLabel) || - (!useMaskLabel && maskImageBufferRange[indexValue] != MaskPixelType{})) && - (confidenceImageBufferRange.empty() || confidenceImageBufferRange[indexValue] > 0.0)) + RealType pixel = std::exp(subtracterImageBufferRange[indexValue]); + N += 1.0; + + if (N > 1.0) { - RealType pixel = std::exp(subtracterImageBufferRange[indexValue]); - N += 1.0; - - if (N > 1.0) - { - sigma = sigma + itk::Math::sqr(pixel - mu) * (N - 1.0) / N; - } - mu = mu * (1.0 - 1.0 / N) + pixel / N; + sigma = sigma + itk::Math::sqr(pixel - mu) * (N - 1.0) / N; } + mu = mu * (1.0 - 1.0 / N) + pixel / N; } - sigma = std::sqrt(sigma / (N - 1.0)); - - return (sigma / mu); } + sigma = std::sqrt(sigma / (N - 1.0)); - template - void N4BiasFieldCorrectionImageFilter::PrintSelf(std::ostream & os, - Indent indent) const - { - Superclass::PrintSelf(os, indent); - - os << indent << "Use Mask Label: " << m_UseMaskLabel << std::endl; - os << indent << "Mask label: " << static_cast::PrintType>(this->m_MaskLabel) - << std::endl; - os << indent << "Number of histogram bins: " << this->m_NumberOfHistogramBins << std::endl; - os << indent << "Wiener filter noise: " << this->m_WienerFilterNoise << std::endl; - os << indent << "Bias field FWHM: " << this->m_BiasFieldFullWidthAtHalfMaximum << std::endl; - os << indent << "Maximum number of iterations: " << this->m_MaximumNumberOfIterations << std::endl; - os << indent << "Convergence threshold: " << this->m_ConvergenceThreshold << std::endl; - os << indent << "Spline order: " << this->m_SplineOrder << std::endl; - os << indent << "Number of fitting levels: " << this->m_NumberOfFittingLevels << std::endl; - os << indent << "Number of control points: " << this->m_NumberOfControlPoints << std::endl; - os << indent << "CurrentConvergenceMeasurement: " << this->m_CurrentConvergenceMeasurement << std::endl; - os << indent << "CurrentLevel: " << this->m_CurrentLevel << std::endl; - os << indent << "ElapsedIterations: " << this->m_ElapsedIterations << std::endl; - itkPrintSelfObjectMacro(LogBiasFieldControlPointLattice); - } + return (sigma / mu); +} + +template +void +N4BiasFieldCorrectionImageFilter::PrintSelf(std::ostream & os, + Indent indent) const +{ + Superclass::PrintSelf(os, indent); + + os << indent << "Use Mask Label: " << m_UseMaskLabel << std::endl; + os << indent << "Mask label: " << static_cast::PrintType>(this->m_MaskLabel) + << std::endl; + os << indent << "Number of histogram bins: " << this->m_NumberOfHistogramBins << std::endl; + os << indent << "Wiener filter noise: " << this->m_WienerFilterNoise << std::endl; + os << indent << "Bias field FWHM: " << this->m_BiasFieldFullWidthAtHalfMaximum << std::endl; + os << indent << "Maximum number of iterations: " << this->m_MaximumNumberOfIterations << std::endl; + os << indent << "Convergence threshold: " << this->m_ConvergenceThreshold << std::endl; + os << indent << "Spline order: " << this->m_SplineOrder << std::endl; + os << indent << "Number of fitting levels: " << this->m_NumberOfFittingLevels << std::endl; + os << indent << "Number of control points: " << this->m_NumberOfControlPoints << std::endl; + os << indent << "CurrentConvergenceMeasurement: " << this->m_CurrentConvergenceMeasurement << std::endl; + os << indent << "CurrentLevel: " << this->m_CurrentLevel << std::endl; + os << indent << "ElapsedIterations: " << this->m_ElapsedIterations << std::endl; + itkPrintSelfObjectMacro(LogBiasFieldControlPointLattice); +} } // end namespace itk diff --git a/Modules/IO/IPL/include/itkIPLFileNameList.h b/Modules/IO/IPL/include/itkIPLFileNameList.h index 94589e5c108c..9b0f3088ba5c 100644 --- a/Modules/IO/IPL/include/itkIPLFileNameList.h +++ b/Modules/IO/IPL/include/itkIPLFileNameList.h @@ -39,11 +39,16 @@ /** Set built-in type. Creates member Set"name"() (e.g., SetVisibility()); */ #define IPLSetMacroDeclaration(name, type) virtual void Set##name(const type _arg); -#define IPLSetMacroDefinition(class, name, type) \ - void class ::Set##name(const type _arg) \ - { \ - CLANG_PRAGMA_PUSH \ - CLANG_SUPPRESS_Wfloat_equal if (this->m_##name != _arg) CLANG_PRAGMA_POP { this->m_##name = _arg; } \ +#define IPLSetMacroDefinition(class, name, type) \ + void class ::Set##name(const type _arg) \ + { \ + GCC_PRAGMA_PUSH \ + GCC_SUPPRESS_Wfloat_equal \ + if (this->m_##name != _arg) \ + { \ + this->m_##name = _arg; \ + } \ + GCC_PRAGMA_POP \ } /** Get built-in type. Creates member Get"name"() (e.g., GetVisibility()); */ diff --git a/Modules/IO/ImageBase/include/itkImageSeriesWriter.hxx b/Modules/IO/ImageBase/include/itkImageSeriesWriter.hxx index 4d82521de87b..e43ae60a608f 100644 --- a/Modules/IO/ImageBase/include/itkImageSeriesWriter.hxx +++ b/Modules/IO/ImageBase/include/itkImageSeriesWriter.hxx @@ -140,7 +140,10 @@ ImageSeriesWriter::GenerateNumericFileNames() for (unsigned int slice = 0; slice < numberOfFiles; ++slice) { + GCC_PRAGMA_PUSH + GCC_SUPPRESS_Wformat_nonliteral snprintf(fileName, IOCommon::ITK_MAXPATHLEN + 1, m_SeriesFormat.c_str(), fileNumber); + GCC_PRAGMA_POP m_FileNames.push_back(fileName); fileNumber += this->m_IncrementIndex; } diff --git a/Modules/IO/ImageBase/src/itkNumericSeriesFileNames.cxx b/Modules/IO/ImageBase/src/itkNumericSeriesFileNames.cxx index d161894b3dca..8fcd873b1d06 100644 --- a/Modules/IO/ImageBase/src/itkNumericSeriesFileNames.cxx +++ b/Modules/IO/ImageBase/src/itkNumericSeriesFileNames.cxx @@ -48,7 +48,10 @@ NumericSeriesFileNames::GetFileNames() // snprintf returns the length of the string generated by the // snprintf operation (excluding the trailing null), // so we can use a 'dry run' to size the buffer. + GCC_PRAGMA_PUSH + GCC_SUPPRESS_Wformat_nonliteral OffsetValueType nchars = snprintf(c, sizeof(c), m_SeriesFormat.c_str(), i); + GCC_PRAGMA_POP if (nchars <= -1) { // broken snprintf reports overflow with -1 @@ -57,7 +60,10 @@ NumericSeriesFileNames::GetFileNames() } OffsetValueType bufflen = nchars + 1; const auto temp = make_unique_for_overwrite(bufflen); + GCC_PRAGMA_PUSH + GCC_SUPPRESS_Wformat_nonliteral OffsetValueType result = snprintf(temp.get(), bufflen, m_SeriesFormat.c_str(), i); + GCC_PRAGMA_POP if (result < 0 || result >= bufflen) { std::stringstream message_cache; diff --git a/Modules/IO/MeshBase/include/itkMeshIOTestHelper.h b/Modules/IO/MeshBase/include/itkMeshIOTestHelper.h index 32c52293d51a..efd6fca48b32 100644 --- a/Modules/IO/MeshBase/include/itkMeshIOTestHelper.h +++ b/Modules/IO/MeshBase/include/itkMeshIOTestHelper.h @@ -39,16 +39,18 @@ // Define a local macro for variable to command testing to avoid including // itkTestingMacros.h, which causes link issues as the module hosting this // file is not a testing module. -#define LOCAL_ITK_TEST_SET_GET_VALUE(variable, command) \ - CLANG_PRAGMA_PUSH \ - CLANG_SUPPRESS_Wfloat_equal if (variable != command) CLANG_PRAGMA_POP \ - { \ - std::cerr << "Error in " << #command << std::endl; \ - std::cerr << " In " __FILE__ ", line " << __LINE__ << std::endl; \ - std::cerr << "Expected " << variable << std::endl; \ - std::cerr << "but got " << command << std::endl; \ - return EXIT_FAILURE; \ - } \ +#define LOCAL_ITK_TEST_SET_GET_VALUE(variable, command) \ + GCC_PRAGMA_PUSH \ + GCC_SUPPRESS_Wfloat_equal \ + if (variable != command) \ + { \ + std::cerr << "Error in " << #command << std::endl; \ + std::cerr << " In " __FILE__ ", line " << __LINE__ << std::endl; \ + std::cerr << "Expected " << variable << std::endl; \ + std::cerr << "but got " << command << std::endl; \ + return EXIT_FAILURE; \ + } \ + GCC_PRAGMA_POP \ ITK_MACROEND_NOOP_STATEMENT diff --git a/Modules/Segmentation/Watersheds/include/itkWatershedImageFilter.hxx b/Modules/Segmentation/Watersheds/include/itkWatershedImageFilter.hxx index 34bcef4fbb58..93498572cfde 100644 --- a/Modules/Segmentation/Watersheds/include/itkWatershedImageFilter.hxx +++ b/Modules/Segmentation/Watersheds/include/itkWatershedImageFilter.hxx @@ -33,8 +33,9 @@ WatershedImageFilter::SetThreshold(double val) val = 1.0; } - CLANG_PRAGMA_PUSH - CLANG_SUPPRESS_Wfloat_equal if (val != m_Threshold) CLANG_PRAGMA_POP + GCC_PRAGMA_PUSH + GCC_SUPPRESS_Wfloat_equal + if (val != m_Threshold) { m_Threshold = val; m_Segmenter->SetThreshold(m_Threshold); @@ -42,6 +43,7 @@ WatershedImageFilter::SetThreshold(double val) m_ThresholdChanged = true; this->Modified(); } + GCC_PRAGMA_POP } template @@ -57,8 +59,9 @@ WatershedImageFilter::SetLevel(double val) val = 1.0; } - CLANG_PRAGMA_PUSH - CLANG_SUPPRESS_Wfloat_equal if (val != m_Level) CLANG_PRAGMA_POP + GCC_PRAGMA_PUSH + GCC_SUPPRESS_Wfloat_equal + if (val != m_Level) { m_Level = val; m_TreeGenerator->SetFloodLevel(m_Level); @@ -67,6 +70,7 @@ WatershedImageFilter::SetLevel(double val) m_LevelChanged = true; this->Modified(); } + GCC_PRAGMA_POP } template diff --git a/Modules/Segmentation/Watersheds/include/itkWatershedSegmenter.hxx b/Modules/Segmentation/Watersheds/include/itkWatershedSegmenter.hxx index fd9e28eeb793..830ca6d0c8f2 100644 --- a/Modules/Segmentation/Watersheds/include/itkWatershedSegmenter.hxx +++ b/Modules/Segmentation/Watersheds/include/itkWatershedSegmenter.hxx @@ -175,18 +175,13 @@ Segmenter::GenerateData() Self::MinMax(input, regionToProcess, minimum, maximum); // cap the maximum in the image so that we can always define a pixel // value that is one greater than the maximum value in the image. - if (std::is_integral_v - // clang-format off -CLANG_PRAGMA_PUSH -CLANG_SUPPRESS_Wfloat_equal - // clang-format on - && maximum == NumericTraits::max()) - // clang-format off -CLANG_PRAGMA_POP - // clang-format on - { - maximum -= NumericTraits::OneValue(); - } + GCC_PRAGMA_PUSH + GCC_SUPPRESS_Wfloat_equal + if (std::is_integral_v && maximum == NumericTraits::max()) + { + maximum -= NumericTraits::OneValue(); + } + GCC_PRAGMA_POP // threshold the image. Self::Threshold(thresholdImage, input, @@ -1201,12 +1196,14 @@ Segmenter::Threshold(InputImageTypePointer destination, { dIt.Set(threshold); } - CLANG_PRAGMA_PUSH - CLANG_SUPPRESS_Wfloat_equal else if (tmp == NumericTraits::max()) CLANG_PRAGMA_POP + GCC_PRAGMA_PUSH + GCC_SUPPRESS_Wfloat_equal + else if (tmp == NumericTraits::max()) { dIt.Set(tmp - NumericTraits::OneValue()); } else { dIt.Set(tmp); } + GCC_PRAGMA_POP ++dIt; ++sIt; }