From a34958585142fe22ec2b0c83409810398944ee55 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Onur=20=C3=9Clgen?= Date: Wed, 28 Jun 2023 18:10:27 +0100 Subject: [PATCH] Fix a bug causing wrong calculation of the affine transformation matrix --- niftyreg_build_version.txt | 2 +- reg-apps/reg_aladin.cpp | 4 ++-- reg-lib/_reg_aladin.cpp | 4 ++-- reg-lib/_reg_aladin_sym.cpp | 5 +---- reg-lib/_reg_base.cpp | 12 ++++++------ reg-lib/cpu/_reg_tools.cpp | 22 +++++++++++----------- 6 files changed, 23 insertions(+), 26 deletions(-) diff --git a/niftyreg_build_version.txt b/niftyreg_build_version.txt index 2b930fc4..c1d1ffbb 100644 --- a/niftyreg_build_version.txt +++ b/niftyreg_build_version.txt @@ -1 +1 @@ -265 +266 diff --git a/reg-apps/reg_aladin.cpp b/reg-apps/reg_aladin.cpp index 1ced15cb..c9c82ec5 100755 --- a/reg-apps/reg_aladin.cpp +++ b/reg-apps/reg_aladin.cpp @@ -163,9 +163,9 @@ int main(int argc, char **argv) { float floatingSigma = 0; float referenceSigma = 0; - float referenceLowerThr = std::numeric_limits::min(); + float referenceLowerThr = std::numeric_limits::lowest(); float referenceUpperThr = std::numeric_limits::max(); - float floatingLowerThr = std::numeric_limits::min(); + float floatingLowerThr = std::numeric_limits::lowest(); float floatingUpperThr = std::numeric_limits::max(); float paddingValue = std::numeric_limits::quiet_NaN(); diff --git a/reg-lib/_reg_aladin.cpp b/reg-lib/_reg_aladin.cpp index 70df10c2..0cc6aa68 100644 --- a/reg-lib/_reg_aladin.cpp +++ b/reg-lib/_reg_aladin.cpp @@ -32,11 +32,11 @@ reg_aladin::reg_aladin() { this->floatingSigma = 0; this->referenceSigma = 0; + this->referenceLowerThreshold = std::numeric_limits::lowest(); this->referenceUpperThreshold = std::numeric_limits::max(); - this->referenceLowerThreshold = std::numeric_limits::min(); + this->floatingLowerThreshold = std::numeric_limits::lowest(); this->floatingUpperThreshold = std::numeric_limits::max(); - this->floatingLowerThreshold = std::numeric_limits::min(); this->warpedPaddingValue = std::numeric_limits::quiet_NaN(); diff --git a/reg-lib/_reg_aladin_sym.cpp b/reg-lib/_reg_aladin_sym.cpp index fe97cca0..f131fea6 100644 --- a/reg-lib/_reg_aladin_sym.cpp +++ b/reg-lib/_reg_aladin_sym.cpp @@ -11,9 +11,6 @@ reg_aladin_sym::reg_aladin_sym() this->backwardBlockMatchingParams = nullptr; - this->floatingUpperThreshold = std::numeric_limits::max(); - this->floatingLowerThreshold = std::numeric_limits::min(); - #ifndef NDEBUG reg_print_msg_debug("reg_aladin_sym constructor called"); #endif @@ -63,7 +60,7 @@ void reg_aladin_sym::InitialiseRegistration() { } } } - if (this->floatingLowerThreshold != std::numeric_limits::min()) { + if (this->floatingLowerThreshold != std::numeric_limits::lowest()) { for (unsigned l = 0; l < this->levelsToPerform; ++l) { T *refPtr = static_cast(this->floatingPyramid[l]->data); int *mskPtr = this->floatingMaskPyramid[l].get(); diff --git a/reg-lib/_reg_base.cpp b/reg-lib/_reg_base.cpp index 8e208d96..54eb63ab 100644 --- a/reg-lib/_reg_base.cpp +++ b/reg-lib/_reg_base.cpp @@ -34,14 +34,14 @@ reg_base::reg_base(int refTimePoint, int floTimePoint) { referenceSmoothingSigma = 0; floatingSmoothingSigma = 0; + referenceThresholdLow.reset(new T[referenceTimePoint]); + std::fill(referenceThresholdLow.get(), referenceThresholdLow.get() + referenceTimePoint, std::numeric_limits::lowest()); referenceThresholdUp.reset(new T[referenceTimePoint]); std::fill(referenceThresholdUp.get(), referenceThresholdUp.get() + referenceTimePoint, std::numeric_limits::max()); - referenceThresholdLow.reset(new T[referenceTimePoint]); - std::fill(referenceThresholdLow.get(), referenceThresholdLow.get() + referenceTimePoint, std::numeric_limits::min()); + floatingThresholdLow.reset(new T[floatingTimePoint]); + std::fill(floatingThresholdLow.get(), floatingThresholdLow.get() + floatingTimePoint, std::numeric_limits::lowest()); floatingThresholdUp.reset(new T[floatingTimePoint]); std::fill(floatingThresholdUp.get(), floatingThresholdUp.get() + floatingTimePoint, std::numeric_limits::max()); - floatingThresholdLow.reset(new T[floatingTimePoint]); - std::fill(floatingThresholdLow.get(), floatingThresholdLow.get() + floatingTimePoint, std::numeric_limits::min()); robustRange = false; warpedPaddingValue = std::numeric_limits::quiet_NaN(); @@ -504,7 +504,7 @@ void reg_base::Initialise() { T *refDataPtr = static_cast(tmpReference->data); reg_heapSort(refDataPtr, tmpReference->nvox); // Update the reference threshold values if no value has been setup by the user - if (referenceThresholdLow[0] == std::numeric_limits::min()) + if (referenceThresholdLow[0] == std::numeric_limits::lowest()) referenceThresholdLow[0] = refDataPtr[(int)reg_round((float)tmpReference->nvox * 0.02f)]; if (referenceThresholdUp[0] == std::numeric_limits::max()) referenceThresholdUp[0] = refDataPtr[(int)reg_round((float)tmpReference->nvox * 0.98f)]; @@ -516,7 +516,7 @@ void reg_base::Initialise() { T *floDataPtr = static_cast(tmpFloating->data); reg_heapSort(floDataPtr, tmpFloating->nvox); // Update the floating threshold values if no value has been setup by the user - if (floatingThresholdLow[0] == std::numeric_limits::min()) + if (floatingThresholdLow[0] == std::numeric_limits::lowest()) floatingThresholdLow[0] = floDataPtr[(int)reg_round((float)tmpFloating->nvox * 0.02f)]; if (floatingThresholdUp[0] == std::numeric_limits::max()) floatingThresholdUp[0] = floDataPtr[(int)reg_round((float)tmpFloating->nvox * 0.98f)]; diff --git a/reg-lib/cpu/_reg_tools.cpp b/reg-lib/cpu/_reg_tools.cpp index 015be4d4..9b4dc6f9 100755 --- a/reg-lib/cpu/_reg_tools.cpp +++ b/reg-lib/cpu/_reg_tools.cpp @@ -104,35 +104,35 @@ void reg_intensityRescale_core(nifti_image *image, switch (image->datatype) { case NIFTI_TYPE_UINT8: currentMin = (DataType)std::numeric_limits::max(); - currentMax = 0; + currentMax = (DataType)std::numeric_limits::lowest(); break; case NIFTI_TYPE_INT8: currentMin = (DataType)std::numeric_limits::max(); - currentMax = (DataType)std::numeric_limits::min(); + currentMax = (DataType)std::numeric_limits::lowest(); break; case NIFTI_TYPE_UINT16: currentMin = (DataType)std::numeric_limits::max(); - currentMax = (DataType)std::numeric_limits::min(); + currentMax = (DataType)std::numeric_limits::lowest(); break; case NIFTI_TYPE_INT16: currentMin = (DataType)std::numeric_limits::max(); - currentMax = (DataType)std::numeric_limits::min(); + currentMax = (DataType)std::numeric_limits::lowest(); break; case NIFTI_TYPE_UINT32: currentMin = (DataType)std::numeric_limits::max(); - currentMax = (DataType)std::numeric_limits::min(); + currentMax = (DataType)std::numeric_limits::lowest(); break; case NIFTI_TYPE_INT32: currentMin = (DataType)std::numeric_limits::max(); - currentMax = (DataType)std::numeric_limits::min(); + currentMax = (DataType)std::numeric_limits::lowest(); break; case NIFTI_TYPE_FLOAT32: currentMin = (DataType)std::numeric_limits::max(); - currentMax = (DataType)std::numeric_limits::min(); + currentMax = (DataType)std::numeric_limits::lowest(); break; case NIFTI_TYPE_FLOAT64: currentMin = (DataType)std::numeric_limits::max(); - currentMax = (DataType)std::numeric_limits::min(); + currentMax = (DataType)std::numeric_limits::lowest(); break; } @@ -284,7 +284,7 @@ template void reg_thresholdImage(nifti_image *image, T lowThr, T upThr) { DataType *imagePtr = static_cast(image->data); T currentMin = std::numeric_limits::max(); - T currentMax = std::numeric_limits::min(); + T currentMax = std::numeric_limits::lowest(); if (image->scl_slope == 0)image->scl_slope = 1.0; @@ -1338,7 +1338,7 @@ void reg_tools_labelKernelConvolution_core(nifti_image *image, } currIterator = tmp_lab.begin(); maxindex = 0; - maxval = std::numeric_limits::min(); + maxval = std::numeric_limits::lowest(); while (currIterator != tmp_lab.end()) { if (currIterator->second > maxval) { maxindex = currIterator->first; @@ -2008,7 +2008,7 @@ DataType reg_tools_getMinMaxValue(const nifti_image *image, int timepoint, bool reg_print_msg_error("reg_tools_getMinMaxValue. The required time point does not exists"); const DataType *imgPtr = static_cast(image->data); - DataType retValue = calcMin ? std::numeric_limits::max() : std::numeric_limits::min(); + DataType retValue = calcMin ? std::numeric_limits::max() : std::numeric_limits::lowest(); const size_t voxelNumber = CalcVoxelNumber(*image); const float sclSlope = image->scl_slope == 0 ? 1 : image->scl_slope;